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ABSTRACT 


As  continued  efforts  are  made  to  reduce  the  radar  cross  sections  of  aircraft  and 
ships,  designs  are  first  modeled  with  computers  and  then  tested  in  the  lab.  In  the  far-field 
of  these  tested  objects,  actual  sources  of  high  reflectivity  or  “Hot  Spots”  on  the  tested 
objects  can  be  isolated  to  within  only  one  half  the  wavelength  of  the  electromagnetic 
wave  used  for  testing.  Ideally,  a  probe  could  measure  fields  on  the  surface  of  the  object 
being  tested  to  completely  isolate  the  source  of  the  hot  spot.  Unfortunately,  the  presence 
of  the  probe  on  the  surface  of  the  object  will  disturb  the  very  fields  it  is  attempting  to 
measure.  Probe  measurements  made  in  the  near  field,  close  to  but  not  on  the  object,  can 
be  designed  to  reduce  the  influence  of  the  probe  while  providing  accurate  field  data.  The 
data  thus  measured,  while  not  able  to  determine  the  source  location  perfectly,  can  be 
used  to  localize  a  source  to  less  than  one  half  wavelength,  the  far-field  diffraction  limit. 
This  thesis  tests  a  technique  for  back  propagating  computer  generated  near  field 
measurements  of  an  axisymmetric  field  source  to  determine  the  fields  closer  to  the 
source.  Several  cases  are  examined  that  test  the  accuracy  and  resolving  capability  of  the 
technique. 
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L  INTRODUCTION 


The  harder  it  is  to  detect  a  platform,  the  easier  it  is  for  a  platform  to  escape 
observation  and  targeting.  This  is  the  driving  force  behind  research  and  engineering  in 
the  area  of  “stealth.”  Stealth  refers  not  only  to  the  oddly  shaped  airframe  of  the  F-1 17, 
but  to  a  wide  realm  of  techniques  designed  to  reduce  the  many  facets  of  a  platform’s 
“signature.”  The  signature  of  a  platform  consists  of  acoustic,  thermal,  and 
electromagnetic  somces  that  are  exploited  in  the  detection  of  that  platform.  As  stealth 
technology  has  become  more  publicized,  various  efforts  for  countering  it  have  also  been 
under  exploration.  Advances  made  in  overcoming  stealth  directly  imply  that  continued 
efforts  are  required  to  further  reduce  signatures. 

A.  GOAL  OF  RESEARCH 

In  the  areas  of  Radar  Cross  Section  (RCS)  and  electromagnetic  sources,  the 
general  trend  has  been  to  attempt  back-propagation  of  field  measurements  from  the  far 
field  region  of  a  platform  down  to  the  surface  of  the  platform  itself  From  this  estimation 
of  the  fields  on  the  platform,  the  distribution,  localization,  and  characterization  of  source 
currents  can  be  approximated.  Unfortunately,  sources  can  only  be  distinguished 
separately  if  they  are  greater  than  one  half  wavelength  apart  when  measurements  are 
taken  from  the  far  field.  Figure  1  shows  the  relationship  of  three  regions  of  interest  to  the 
surface  of  a  source  and  the  regions’  resolution  limitations. 

Alternatively,  since  exact  source  localization  would  be  ideal,  field  measurements 
made  directly  on  the  platform  by  a  perfect  probe  might  yield  the  required  information. 
Unfortunately,  the  presence  of  an  actual  probe  will  greatly  disturb  the  fields  and  source 
distribution  on  the  platform. 

This  leaves  the  near  field  region  of  the  platform,  between  the  platform’s  surface 
and  the  far  field,  in  which  to  make  measurements  and  to  image  the  sources.  For  a  probe 
in  the  near  field,  the  measurable  fields  no  longer  constrain  the  source  resolution  to 
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structures  greater  than  one  half  wavelength  (as  was  the  case  for  the  far  field)  and  allow 
for  greater  refinement.  Additionally,  the  probe  can  be  placed  at  a  distance  sufficiently  far 
away  from  the  surface  to  make  the  effects  of  its  presence  managable. 


Figure  1.  Far  Field/Near  Field/Surface  Localization  Limitations 

Exploring  a  technique  for  the  back  propagation  of  fields  measmed  in  the  near 
field  region  at  one  distance  from  the  source  onto  fields  closer  to  the  source  is  the  primary 
goal  of  this  thesis. 

B.  BACKGROUND  OF  RESEARCH 

The  analysis  of  near  field  phenomena  that  spurred  the  current  research  originated 
in  the  field  of  acoustics.  For  plate  and  cylinder  like  structures,  a  newly  defined  quantity 
called  the  supersonic  acoustic  intensity  vector  provided  a  tool  for  the  location  and 
characterization  of  acoustic  sources  from  near  field  acoustic  measurements.  [Ref  1] 
From  acoustics,  the  supersonic  analysis  of  sound  wave  components  was  translated  into 
the  realm  of  electromagnetics.  For  cylindrical  geometiy,  the  now  superluminal  modes 
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which  satisfy  <  A:  and  have  faster-than-light  propagation  along  the  cylindrical  axis 

provide  time-average  power  to  the  far  field.  Measurements  taken  in  the  near  field  were 
back-propagated  using  an  optimal  deconvolution  filter  to  account  for  the  subluminal 
modes.  [Ref  2]  This  technique  of  using  superluminal  modes  was  then  extended  fi'om 
separable  geometries  to  non-separable  geometries.  Using  the  Finite  Element  Method  to 
relate  the  boimdaiy  conditions  of  a  measured  field  to  the  surface  currents  induced  on  an 
axisymmetric  body  by  scattering,  fields  were  back-propagated  to  estimate  the  source 
distribution.  Resonance  issues  in  the  spherical  test  case  complicated  the  overall 
effectiveness  of  that  particular  approach.  [Ref  3] 

This  thesis  continues  to  explore  methods  of  near  field  back-propagation  that  apply 
to  general  non-separable  geometries  from  axisymmetric  bodies.  These  are  bodies  that  do 
not  have  explicit  separable  modes  such  as  cylindrical  or  spherical  wave  functions. 
Unlike  the  previous  work  in  non-separable  geometries,  this  thesis  attempts  to  reconstruct 
fields  closer  to  the  source  and  not  strictly  the  source  itself  Additionally,  this  thesis  looks 
at  source  current  distributions  that  are  not  induced  by  scattering. 

This  study  is  organized  as  follows:  Various  cases  of  generated  electromagnetic 
radiation  source  fields  are  discussed  in  Chapter  H.  Various  techniques  in  calculating 
forward  propagation  fields  and  testing  their  convergence  are  examined  in  Chapter  HI. 
Various  techniques  in  calculating  back  propagation  of  fields  are  discussed  in  Chapter  IV. 
A  summary  and  concluding  discussion  are  presented  in  Chapter  V.  An  appendix  has 
been  included  with  the  fundamental  set  of  computer  codes  to  allow  for  an  indepth  review 
and  continuation  of  this  line  of  research. 
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n.  GENERATION  OF  ELECTROMAGNETIC  SOURCE  FIELDS  AND  THE 


FIELD  TRANSFORMATION  OPERATOR 

A.  GENERATION  OF  THE  EXACT  FIELDS 


Before  one  can  perform  numerical  experiments  on  the  accuracy  of  back- 
propagated  fields,  one  must  generate  the  fields  to  be  propagated  and  transform  matrix 
that  does  the  propagation  operations.  First  we  consider  a  single  dipole  located  at  the 
origin.  Following  the  derivation  [Ref  4],  a  closed  form  solution  to  the  radiation  patterns 
any  distance  away  from  the  source  exists.  These  closed  form  solutions  form  the  exact 
fields  at  various  p  values  which  are  in  turn  used  for  two  purposes.  First,  the  exact  fields 
are  used  as  checks  of  the  accuracy  for  various  numerical  techniques,  and  second,  as  the 
basis  for  determining  the  eqmvalent  electric  and  magnetic  currents  required  for 
numerical  propagation.  The  dipole  is  placed  at  the  origin  with  its  length  along  the  z-axis 
as  shown  in  Figure  2. 
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The  axisymmetric  fields  it  produces  will  be  the  same  at  any  point  with  the  same  p 
and  z  around  the  axis.  Any  combination  of  dipoles  oriented  along  the  z-axis  will  also 
produce  axisymmetric  fields.  The  closed  form  solutions  are  based  on  an  analysis  of  a 
dipole  with  a  sinusoidal  current  distribution.  Figure  3  shows  the  dipole  geometry  used  in 
the  derivation.  The  current  distribution  on  the  dipole  is  sinusoidual  such  that 


Figure  3.  Dipole  Exact  Field  Derivation  Geometry 

Ig  =  where  z'  €  \yh-h\.  Figure  4  illustrates  two  cases.  The  first 

case  shows  the  dipole  with  an  overall  length  less  then  one  half  the  wavelength  associated 
with  the  current  fi-equency.  The  second  case  shows  the  dipole  with  a  length  between  one 
half  and  one  wavelength.  Phaser  Fields  may  be  determined  through  use  of  the  vector 
potential,  computed  for  line  currents  as  [Ref  5] 

~~jkR 

A{x,y,z)  =  -^jlXx\y,z’)^-—dz’  (1) 

with  R  =  +y^  +{z-  /f  .  (2) 
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Current  Distrubutions  for  Several  Frequencies  with  Dipole  Length  <=  X/2;  Dipole  Length  =  0.1 


Current  Distrubutions  for  Several  Frequencies  with  Dipole  Length  >  )J2\  Dipole  Length  =  0.1 


y  in  meters  for  L  =  0.1 


Figure  4.  Sinusodial  Current  Distribution 
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Shifting  to  cylindrical  coordinates  and  substituting  the  sinusoidal  current  distribution  of 
the  dipole  for  the  integration  may  be  readily  performed.  Using  Maxwell’s 

equations,  the  magnetic  field  due  to  vector  potential  A{p,^,z^  is  given  by 


H^=—^xa).  The  resulting  electric  field  is  given  by  £"  =  — — (v x ir).  The  fields 
produced  by  a  sinusoidal  current  distribution  can  then  be  derived  as  the  following, 


H.  =  _(2cos)t/z)e-"^] 

jAnp 


~  DIPOLE 


Jz-h] 

e~MA 

-J- 

(z  +  h\ 

A  P  J 

1  P  J 

R. 


■(2cosM) 


z  e 


■M, 


P  K 


Rb 


R. 


-(2cosM)- 


R. 


(3) 

(4) 

(5) 


These  three  equations  are  functions  of 

1)  Dipole  height,  h, 

2)  Wavelength,  k  =  , 

c 

3)  Maximum  Current  Magnitude,  Idipole,  and  the 

4)  Observation  point  in  p  and  z,  but  not  (j)  since  the  fields  are  axisymmetric. 
Using  H<j,  and  E^,  the  fields  tangential  to  a  cylindrical  surface  of  any  radius  p  can  be 
determined  within  machine  roimd  off  accuracy. 


B.  TANGENTIAL  FIELD  GENERATION  FROM  EXACT  AXISYMMETRIC 
FIELDS  ON  A  CYLINDRICAL  SURFACE 


The  fields  on  a  cylindrical  surface  at  pz  can  be  determined  in  two  ways.  First,  the 
fields  can  be  calculated  fi'om  the  actual  source  axisymmetric  current  distribution. 
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Second,  if  the  fields  are  exactly  known  on  some  closed  surface  that  contains  the  actual 
sources  and  is  inside  the  p2  cylindrical  surface,  then  the  surface  equivalence  theorem  can 
be  used  to  transform  the  known  fields  into  electric  and  magnetic  currents  on  the  closed 
surface.  Since  the  source  current  distribution  is  axisymmetric,  the  inner  surface  is  chosen 
to  be  an  infinite  cylindrical  surface.  Figure  5  shows  the  relative  geometry  of  the  two 
coaxial  cylindrical  surfaces  used  in  determining  the  fields  at  p2. 


Figure  5.  Coaxial  Cylindrical  Surfaces  Geometry 

The  exact  fields  from  the  dipole  source  have  both  E-  and  H-field  components 
tangential  to  a  cylindrical  surface  as  shown  in  equations  (3)  and  (5).  Applying  the 
surface  equivalence  theorem,  the  fields  at  pi^  become  the  following  equivalent  surface 
currents. 

=^xS  =  =  =JA  («) 
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(7) 


Using  J surface  and  M surface,  thc  vectof  potentials  A  and  F  can  be  determined  via  the 
inhomogeneous  vector  potential  equations. 


A  = 


A 


- 

- - ds’ 

4^)jj  /? 


F  =  ^((M^^,^-^-ds' 
4^JJ  o 


jtR 


R 


(8) 

(9) 


Substituting  A  and  F  into  Maxwell’s  equations  yields  the  following  component  fields 
[Ref5]. 


Ha  = 


Ea  = 


( 


1 


WxHa 


1 


Ef  = - V  X  F 


Hf  = - —VxEf 


H  =  Ha^Hf 
E  =  Ea  +  Ef 


(10) 

(11) 

(12) 

(13) 

(14) 


—  — 

In  order  to  determine  if  and  E,  ultimately  the  Green’s  kernel  -  will  have  to  be 

R 

integrated.  The  following  subsection  sets  up  the  problem  for  that  integration. 
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1.  Setting  Up  the  Problem 


The  geometry  for  an  arbitrary  field  point  is  shown  in  Figure  6. 


Figure  6.  General  Case  for  an  Arbitrary  Field  Point  (pz,  ^2,  Za) 

The  steps  for  setting  up  the  problem  for  an  arbitrary  field  point  are  shown  below. 


p  =  cos  ^2^  +  sin  ^2y 

(15) 

^  =  cos(90  +  ij)^  +  sin(90  + 

(16) 

p'  =  cos^x  +  sin^y 

(17) 

=  cos(90  +  ^j)x  +  sin(90  +  <f)^y 

(18) 
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A  A# 

z-z 

(19) 

r  —  p^p  + 

(20) 

r^p^p  +  z^z 

(21) 

R  =  r-r'  =  p^p- p^p’  +  (zj  -  Zi)z 

(22) 

R  =  \p2- Pi - <^2)1^  +  \pi{- sin(A - +  (^2 “ 

(23) 

J  =  J^z'  =  J^z 

(24) 

M  =  M/  =  sin(^J,  -  ^^2))]^  +  [m^(cos((1\  -  ^>2))^ 

(25) 

=  ^pI  +  p\  -  Ip^p^  cos(^^,  -  ^^2)+  fe  - 

(26) 

There  does  not  appear  to  be  a  simple  closed  form  solution  to  the  Green’s  kernel 

integration  over  the  infinite  cylindrical  surface  required  to  find  A  and  F .  However, 
there  is  a  straightforward  numerical  technique  for  solving  the  integration.  First  one 
divides  the  cylinder  at  pi  into  many  small  bands  or  rings  that  have  a  constant  current 
magnitude  on  the  ring  as  shown  in  Figure  7.  Then  the  ring  is  subdivided  into  many 
smaller  squares  such  that  both  electric  and  magnetic  equivalent  currents  are  constant  on 
the  square  [Figure  8]. 
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Az  corresponds  to  the  width 
of  the  thin  ring  such  that 

ISz  ^  ISz 

Z. - <  z.  <  Z-  H - 

'  2  ‘  '  2 


Figure  7.  Physical  Meaning  of  Az 


Figure  8.  Subdivision  of  the  Ring  Showing  Constant  Equivalent  Currents 


The  influence  of  the  currents  in  each  square  on  a  ring  to  the  field  at  pz  will  be 
simimed  together  to  give  the  total  influence  of  the  ring.  Since  the  fields  at  all  p  are 
axially  symmetric,  the  influence  of  whole  rings  can  be  used  when  setting  up  the 
mathematical  (matrix)  representation  of  the  system,  as  opposed  to  the  individual 
influences  of  each  small  square  on  each  ring  separately.  One  sums  the  respective 
contributions  to  an  individual  point  on  the  pz  cylindrical  surface  made  by  each  ring  along 
the  length  of  the  iimer  cylinder  to  get  an  estimate  of  the  actual  field  at  that  point  on  pz. 
This  is  repeated  for  each  point  on  pz  in  order  to  obtain  the  estimated  fields  on  the  entire 
Pz  cylinder. 

Taking  advantage  of  the  axial  symmetry  of  the  fields,  the  field  point  to  be 
estimated  is  placed  on  the  x-z  plane,  where  y=0  and  (|)2=0  [Figure  9]. [Ref  6]  When  the 
field  point  is  on  the  x-z  plane,  the  problem  set  up  is  simplified. 


A  A 

p  =  x 

(27) 

II 

(28) 

f  =  Z 

(29) 

M  =  (m^  sin(-  ))yO + (m^  cos  <p)^ 

(30) 

^  =  1/^2  -  A  cos^ilo  +  [/>,(-  sin^i, +  [z^  -  ZjJz 

(31) 

R  =  yjpl  +  pf  -  2p^p2 COS(4  +  (z2  -  zj 

(32) 
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Figure  9.  Field  Point  (xa,  0,  Z2)  and  the  Cylinder  Segment 


By  using  a  ring  of  width  Azi,  A  is  calculated  by  the  summation  of  AA2  segments. 


4^  =  — 


EJW, 


z- - p^d^^hz 


Ri 


(33) 


Therefore,  the  i*^  segment  of  Az=(AA2)i  is  as  follows. 


(34) 


00  ^ 

In  the  limit  of  Az  ->  dz,  AAz  ->  dAz,  and  ->  jdA^  =  A^. 


AA, 


fi.PiAz 

Arc 


"  7/CK 


(35) 


Following  a  similar  set  up  F  and  AF  are  found  to  be  the  following. 

^  =  i  —P,d4,,Az  (36) 

i=-oo_- 
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AF  =  I M,  (37) 

An  R 

Once  the  vector  potentials  have  been  broken  into  segments  for  integration,  the 
determination  for  the  electric  and  magnetic  field  segments  to  be  summed  over  the  inner 
and  outer  cylinders  is  the  next  set  of  steps. 


2.  The  Equivalent  Electric  Current  Case  and  the  Vector  A  Potential 


Using  as  found  in  equation  (35),  can  be  determined  by  the  following 
equation. 

AH7  =  — VxX4  =  — (38) 

ft  ft  I  4;r  {'■  'r  ”J 

Since  the  curl  operation  is  with  regard  to  the  coordinate  system  of  the  observer  (p2,  <t>2,  Z2) 
and  the  integration  is  with  respect  to  (t)i,  then  the  curl  can  be  brought  inside  the 
integration  as  shown  below. 

4^  4L  V  ^  JJ 


p.Az’l 


■.J,i  +  _(Vx/,2)  U, 


Since  JJ:  is  a  function  dependent  only  on  pi  and  Zi  and  the  curl  operation  is  with  regard 
to  the  observer’s  coordinate  system  (p2,  ^2),  then  Y  xj^z  =  0 .  One  will  encounter 

many  varieties  of  the  operation  shown  below.  [Ref  6] 

V—  =  ( (40) 


This  current  case 


-jkR  _  .r  \  _ 

yields  the  following  for  n=  1 ,  V - =  ~  +  “a  e  ■’^R. 

R  \R  R  j 
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Therefore,  one  can  say  the  following. 

(AHa\  =  ^^\P2  +  -  A  jfy +  ]  (44) 

— ;r\  y  y 

{^X  =  A  ]  (45) 

Since  R  =  ■yjpl  +  pi  -lp^P2 cos^  +  (zj  - 2.^  ,  it  is  an  even  function  in  (j)i  over 
the  range  -71  to  +71.  All  functions  of  only  R  are  even.  All  functions  of  R  and  cos(j)i  are 
also  even.  But,  since  sin(j)i  is  an  odd  function,  even  functions  of  R  and  odd  functions  of 
sin(j)i  will  be  odd  over  the  interval  and  will  integrate  to  zero.  Therefore,  (AHa)p  =  0. 
Because  one  will  encounter  various  versions  of  the  functions  of  R  and  combinations  of 
functions  of  R  and  cos(j)i  that  are  to  be  integrated,  one  defines  Un  and  as  well  as  Wn 
(which  will  be  used  in  the  next  section),  as  follows.  [Ref.  6] 


-"Tt 

(46) 

V„=-  jcos^—d^ 

(47) 

W„=l(  cos^J, 

(48) 

Un,  Vn,  and  Wn  are  solved  numerically  by  dividing  the  ring  into  many  segments  [Figure 
5]  and  summing  their  individual  components.  The  Matlab  programs  INTPHI*.M  does 
this  numerical  integration  step.  Since  Jz  and  M<|,  are  uniform  aroimd  the  ring  and  the 
differences  arise  from  the  physical  distance  to  each  segment  and  cos(l)i,  then  the  segments 
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can  be  summed  around  the  ring.  The  ring  will  be  the  basic  element  for  determining 
fields  as  opposed  to  each  A(j)ixAz  segment.  Overall,  this  reduces  the  computation  time 
and  simplifies  the  problem  set  up. 

Substituting  Un,  V„,  and  with  applicable  values  for  n  yields  the  follovwng. 

(t^;),=^^\pAU,+JkU,)-p,(V,  +  JkV^)]  (49) 


3.  The  Equivalent  Magnetic  Current  Case  and  Vector  F  Potential 


The  mathematical  details  of  the  equivalent  magnetic  current  case  for  Ha  and  the 
equivalent  electric  current  for  Ep  are  very  similar  in  form. 


Aff  VxAF 


(50) 


In  AA ,  there  was  only  a  f  component  in  that  contributed  to  AH^ .  In  AF  there  are 
two  components  of  :  a  p  component  and  a  ^  component.  These  two  components  add 


some  complexity  to  the  previous  mathematical  derivation  of  AH^ . 
AAz  f 


Ak 


1 


Vx 


\  \ 


R 


JJ 


d<j>^ 


tiEp  = 


4^ 


I 


R 


xM, 


dA 


p^lSz 


AF_  =  f 

4n  hR^  R^  )  R 


^xM^)dA 


(51) 


(52) 


(53) 


Looking  at  the  Rx  term  gives  the  following  result 

RxM^=m)^{z^-  ZjXcos^,)^  -  (^2  -  ZiXsin^i)^  +  {p^  co&Ai  -  Pi)z\  (54) 


Substituting  /?xM,  into  AE^  and  breaking  it  into  components  yields  the  following 
result  [Ref  6] 
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(55) 


(56) 


(57) 

(58) 

(59) 

4.  Cross  Field  Influence  Term 

In  Section  n.B.2  “The  Equivalent  Electric  Current  Case  and  the  Vector  A 
Potential”  set  fourth  above,  AH^  at  p2  was  calculated  from  AA,  which  was  calculated 
from  J  which,  in  turn,  was  determined  from  a  knowledge  of  the  H-field  at  pi. 
Therefore,  AH^  is  the  part  of  the  total  AH  at  p2  determined  strictly  from  the  H-field  at 
pi-  Similarly  in  Section  in.B.3  “The  Equivalent  Magnetic  Current  Case  and  Vector 
F  Potential,”  AEp-  at  p2  was  calculated  from  AF  ,  which  was  calculated  from  M , 

which  was  determined  from  knowledge  of  the  E-field  at  pi.  Therefore,  AJEp-  is  the  part 
of  the  total  AF  at  P2  determined  strictly  from  the  E-field  at  pi.  But  these  are  only  parts 
of  their  respective  total  fields.  Since  AH  ^otm.  =  ^ a  ^  total  =  +  AF^  , 

the  influence  on  AH  from  the  E-field  at  pi  and  the  influence  on  AE  from  the  H-field  at  pi 


M^p^Az 


(af^)^  - 


_  M^p^Az 


An 


^2  ~ 


(A£,X  =  a(c/3  +  m)] 


(A£A  = 


4;r 
M^p^Az 


An 


=0 
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must  be  determined  in  order  to  accurately  predict  the  fields  at  p2.  Three  ways  to 
determine  these  field  components  were  calculated. 


a.  Method  A:  Maxwell’s  Equations  Applied  Directly  to  Calculated 
AH^  and  AE^  Fields 

This  method,  the  most  straightforward,  takes  the  three  previously 
calculated  equations  of  AH^  and  AEp  and  computes  A£^  and  AH^  from  them.  Using 
Maxwell’s  equations  and  substituting  the  calculated  AH^ ,  one  can  state  the  following. 


1 


AE.  =-^VxAH  =  — 5— VxAH> 
jas^  j(oe^ 


AE,= 


1  r  aAH, .  1  a 


p+- 


ja>s,  [  azj 


(aeA  = 


jcoeM 


\P2 


dU,  dU. 


+  jk- 


\  5^2  5^2  J 


ydz^  dzjj 


(AEj,  =  -^^{2[u,  +  y^F2]+ 

joasATt  p2 


Pi 


ik^ 

ap2  5^2 


-A 


^Pl  ^Pl 


(60) 

(61) 


(62) 


(63) 


Since  all  the  representations  contain  partial  derivatives  of  Un  and  Vq  with  respect  to  p2 
and  Z2,  it  is  useful  to  calculate  them  in  terms  of  other  Un,  Vn,  and  Wn  functions  as  shown 
below.  [Ref  6] 

+ yw„,)  (64) 

^Pi 

-^V,  =  pknW„^^  jkW„y  jkV„,^  (65) 

dp^ 
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(66) 


=  -fe  +y^f4+i) 

&2 

&2 

Substituting  the  partial  derivatives  of  Un  and  Vn  into  AE^ 
following  results. 

+  3yW,  -  eu,\ 

^  j(os,An 

(aeJ.  =  + yW2]-^[E3 + jkv,] 

JCOSM  P2 

+  pl]lc^U,-ijkU,-iU,\ 

+  P  ,P^^V,+6jkV,-lk%\ 

For  the  AH^  fields  the  following  can  be  stated. 

—  1  r7  1  5AE 

AH^  = - V  X  AE^  = - \ - - - 

76)//,  j(Ofi,  [  dz^  dp2 

^  jap^An  {dz^ 

+-^\p2(y,+ jtv,)- p,(u,*  jku,)i 
SPt  J 


(67) 

and  simplifying  yields  the 

(68) 


(69) 

(70) 

(71) 
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Substituting  the  partial  derivatives  of  Un  and  Vn  into  AH^  and  simplifying  yields  the 
following.  [Ref  6] 


^(^3  +/^2  +(^2  +31^^) 

^  jcouM 

+  A a(3^5  +  3 jkW^  - k%  +3Us+3 jkU^  - k^u)  }  (72) 


Two  additional  methods  (Method  B-Maxwell’s  curl  equations  applied  twice  to 
original  vector  potentials  and  Method  C--solving  AE^  and  AH  directly  from  vector 
potentials  by  an  alternative  method)  for  determining  the  cross-field  terms  are  derived  in 
Appendix-A. 
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in.  CALCULATING  FORWARD  PROPAGATION  FIELDS  AND  TESTING 

THEIR  CONVERGENCE 


In  this  chapter,  the  AE  and  AH  segments  are  used  to  generate  a  transfer  matrix. 
The  fields  on  a  cylindrical  surface  at  a  distance  of  p2  from  the  source  can  be  determined 
by  applying  this  transfer  matrix  to  the  equivalent  currents  generated  by  fields  at  a 
distance  of  pi  from  the  source.  The  convergence  of  the  three  methods  for  determining 
AE  and  AH  is  examined.  Choosing  Method  A,  a  parametric  analysis  of  the  convergence 
of  fields  at  P2  is  conducted.  Trends  of  the  convergence  of  the  calculated  fields  to  the 
exact  fields  at  P2  are  examined  based  upon  differing  input  parameters  in  Method  A. 


A.  CALCULATION  OF  TRANSFER  MATRIX  FOR  CALCULATING  THE  pz 
FIELDS 

Each  E^,  Ep ,  zxAHp  segment  depends  on  several  parameters,  the  most 
important  of  which  is  the  distance  R  between  the  observation  point  on  a  cylindrical 
surface  at  a  distance  P2  from  the  source  and  points  on  a  cylindrical  surface  at  pi.  In  order 
to  estimate  the  fields  at  the  i*  point  on  P2,  the  contributions  from  all  points  on  pi  must  be 
summed.  Figure  10  shows  the  physical  relationship  between  the  points  on  pa  and  the 
points  on  pj. 
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Figure  10.  The  Physical  Relationship  Between  Points  on  P2  and  pi 


This  implies  a  vector  multiplication  structure  (as  follows). 


T 

T 


\pl^p2 


(73) 


For  multiple  points  on  P2 ,  the  same  vector  of  fields  at  points  on  pi  is  used  but  each  point 
on  pi  will  have  an  individual  transform  vector  as  shown  below. 


P\->P2 


'Pi-^Pi 


(74) 


(75) 
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If  the  point  distributions  on  pi  and  p2  are  imiform,  then  all  and  all 

.  If  the  point  distributions  are  uniform  and  equivalent  for  both  pi  and  pz, 
then  Az^  =  Az„  =  Az .  If  either  of  the  distributions  on  the  pi  or  Oj  surfaces  is  not 

uniformly  spaced,  then  the  resultant  transfer  matrix  T  will  have  no  commonly 
identifiable  structure.  If  the  point  distributions  are  both  uniformly  spaced,  but  with 
Az„  56  Az„  ,  then  T  will  have  some  structure.  When  Az„  =  Az^  ,  then  the  transfer 

matrix  will  be  banded  and  have  a  Toeplitz  structure  [Figure  1  l][Ref  7],  with  r„  ^  t_„ .  If 
Az„  =  Az^  and  z„  =  ,  then  the  transfer  matrix  will  have  a  Toeplitz  structure  with 

Pi  P\  Pi  Pi  ^  ^ 


Figure  11.  General  Structure  of  a  Toeplitz  Marix 
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The  bands  in  the  transfer  matrix  are  caused  by  the  ^  terms  in  R,  (AE)^ 

and  AH,  and  (z^^^  -  z^^„)  ^  terms  in  (AE)p.  When  Az’s  are  uniform  and  equal,  but  the 

actual  z  distributions  are  not  equivalent  (i.e.  offset  from  each  other),  then  the  transfer 
matrix  elements  are  translation  invariant  and  depend  only  on  m-n.  When  Az’s  are 
uniform  and  the  z  distributions  are  aligned,  then  only  the  (AE)^  term  is  translationally 

invariant  depending  on  m-n.  All  the  other  terms  are  based  on  -  z^  „)  so  they  are 

translationally  invariant  depending  on  1  m-n  | .  Transfer  matrices  with  the  m-n 
dependence  are  full,  nonsymmetric  Toeplitz  matrices,  in  Figure  11.  Transfer 

matrices  with  the  |m-ni  dependence  are  full,  symmetric,  non-Hermetian  Toeplitz 
matrices,  t„  in  Figure  11.  An  interesting  side  note  is  that  transfer  matrices  for 

Azi=CAz2  or  AziI>=Az2,  with  C  and  D  integer  constants  greater  than  one,  are  made  by 
selecting  specific  rows  or  columns  of  a  Toeplitz  matrix  based  on  the  Az  with  the  smallest 
size. 


The  specific  fields  required  for  this  thesis’  study  of  back  propagation  are  the 
and  H<j,  fields  on  the  p2  cylindrical  surface.  These  total  fields  are  foimd  using  a  large 
(2Mx2N)  transfer  matrix  consisting  of  four  blocks.  Each  block  is  the  (MxN)  transfer 
matrix  for  a  specific  contribution  to  the  total  field  from  one  of  the  two  equivalent 
currents  on  the  pi  cylindrical  surface. 


T. 


2Mx2N 


T  T 

T  T 
/he  ^hh. 


[e.h,It  =  [e,h,1 


(79) 

(80) 


In  general  the  complex  elements  of  The  ate  very  large  in  magnitude  compared  to  the  rest 
of  the  matrix  elements.  Conversely,  the  elements  of  Teh  are  very  small  compared  to  the 
other  elements.  The  relative  magnitude  of  the  Tee  and  Thh  elements  are  roughly  the 
same.  This  difference  in  relative  magnitudes  is  due  the  small  H-field  at  pi  contributing  to 
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the  large  E-field  at  p2  for  The,  and,  conversely,  the  large  E-field  at  pi  contributing  to  the 
small  H-field  at  pz  for  Teh. 

Additionally,  T  will  generally  be  rectangular  since  more  points  will  be 
determined  at  pz  than  are  measured  at  pi.  This  is  in  preparation  for  over  determined 
inverse  systems. 


B.  COMPARISON  OF  THREE  METHODS  FOR  DETERMINING  THE 

CALCULATED  FIELDS  AT  pz 

The  three  different  methods  for  determining  the  calculated  fields  at  pz  are 
compared  using  the  Least  Squares  Error  (LSE)  of  each  method’s  field  compared  to  the 
exact  field. 

ISE,  =  (81) 

I  EXACT / 

The  LSE  is  a  measure  of  a  particular  method’s  total  energy  contained  in  the  error  with 
respect  to  the  total  energy  contained  in  the  exact  field,  expressed  as  a  percentage.  A 
comparison  of  the  three  methods  LSE  shown  in  Figure  12  reveals  that  their  differences 
were  indistinguishable.  This  indicates  that  no  particular  method  is  better  than  the  others 
in  the  near  field.  Method  A  was  chosen  for  all  further  investigations.  For  the  fields 
calculated  in  Figure  12  and  in  all  subsequent  figures,  the  frequency  of  the  current 
distribution  on  the  dipole  is  f-  300  MHz.  This  frequency  corresponds  to  a  wavelength 
X,  =  1  meter.  Thus,  values  of  p  and  z  also  represent  values  relative  to  X. 
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Figure  12.  LSEs  for  Three  Methods  on  the  Same  Plot 
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C.  CUMULATIVE  TOTAL  ENERGY  DENSITY  OF  FIELDS 


In  order  to  get  an  idea  about  the  spread  of  energy  contained  in  the  fields  of  a 
single  dipole,  the  cumulative  total  energy  distribution  for  both  Ez  and  H<j,  is  determined  in 
Figure  13  and  Figure  14.  For  a  dipole  centered  at  the  origin^  the  plots  illustrate  the 
percent  of  total  field  energy  in  the  field  for  |  z  |  >Zc ,  the  z  truncation.  These  plots  provide 
a  guideline  for  choosing  a  maximum  z  in  order  to  keep  a  standard  amoirnt  of  field  energy 
when  doing  calculations.  Additionally,  it  is  also  a  measure  of  the  rapidness  of  change  in 
percent  of  energy  as  Zc  approaches  zero.  These  two  considerations  provide  a  starting 
point  for  understanding  the  effects  of  changing  parameters  to  the  amoimt  of  error  to  be 
expected  in  the  calculated  fields  at  pa. 

D.  PARAMETRIC  ANALYSIS  OF  THE  CALCULATED  FIELDS  USING 
METHOD  A 

The  following  subsections  outline  trends  in  the  accuracy  of  the  calculated  fields 
with  respect  to  the  exact  fields  on  the  cylindrical  surface  at  p2.  All  calculations  use 
uniformly  distributed,  equivalent  point  distributions  for  and  z^j » 

maximum/minimum  equal  to  the  z^^  maximum/minimum.  The  relative  merit  of  various 

parameter  changes  can  be  reviewed  using  three  techniques.  First,  the  Real  and  Imaginary 
components  of  the  exact  and  calculated  fields  are  plotted  together.  Second,  the  Relative 
Error  of  the  magnitude  of  the  calculated  field  with  respect  to  the  magnitude  of  the  exact 
field  on  a  point  to  point  basis  is  plotted.  Last,  the  Least  Squares  Error  expressed  as  a 
percentage  is  plotted. 
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Contour  Plots  of  Cumlative  Energy  Distribution  for  Mag(Ep  from  Dipole  at  the  Origin 


Figure  13.  E  Plot 
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p  (m) 


Cutoff-  z  (m);  zMax  =  20  (m);  Dipole  Length  =  0.1  (m);  Nz  =  20000 


Figure  14.  H  Plot 


1.  Varying  the  Truncation  in  Z 


To  get  an  appreciation  of  the  effect  of  var5dng  the  truncation  of  z  at  a  fixed  p 
value,  Figure  15  shows  a  plot  of  the  LSE  in  the  calculated  p2  fields  for  a  range  of 
maximum  z.  From  Figures  13  and  14  it  is  clear  that  as  the  cylindrical  surface  pi  moves 
further  from  the  source  the  value  of  the  z  cutoff  point  must  be  increased  to  maintain  a 
constant  level  of  the  field’s  energy  for  pi  calculations.  Figure  13  shows  that  changes  in 
the  z  truncation  point  have  a  significant  impact  on  the  amount  of  error  in  the  p2  fields 
when  the  magnitude  of  the  maximum  z  is  small.  Then,  beyond  a  certain  point,  increasing 
the  truncation  point  in  z  does  little  to  decrease  the  error  in  magnitude  z  p2  fields.  Figure 
15  confirms  these  two  observations.  For  z  cutoff  less  then  two  meters,  the  LSE  is  greater 
than  one  percent. 

2.  Varying  the  Length  of  Az 

In  building  the  transfer  matrix  and  defining  the  fields  on  the  two  cylindrical 
surfaces,  an  estimation  of  the  actual  smooth  continuous  fields  was  made.  The  fields  were 
estimated  to  be  the  same  over  a  length  of  z,  called  Az.  As  Az  decreases,  higher  spatial 
frequencies  for  fields  can  be  resolved  uniquely.  In  other  words,  the  estimate  of  the 
continuous  field  made  by  discrete  steps  Az  apart  becomes  more  accurate  as  Az  decreases. 
If  there  are  no  spatial  frequencies  higher  than  some  kz,  then  decreasing  Az  when 

'2.71 

Az  «  —  will  provide  little  improvement  to  the  calculated  field  accuracy.  Conversely, 

k^c 
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decreasing  Az  when  Az  »  —  will  dramatically  improve  the  field  accuracy. 

k/: 

These  trends  are  easily  seen  for  cases  when  the  z  truncation  point  is  large  enough 
to  include  most  of  the  energy  in  the  field.  As  seen  in  Figure  16,  when  coarse  estimates  of 


32 


the  field  are  refined,  large  Az  is  decreased  and  the  error  in  the  calculated  fields  also 
decreases. 
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Least  Squares  Error  (%)  for  Mag(Ej)  for  Varying  Max  2  and  fj2 ;  pi  =  0.5  m:  100 


Max  z  (m) :  Nz/meter  =  32 


Least  Squares  Error  (%)  for  Mag(H^  for  Varying  Max  2  and  p2 ;  pi  =  0.5  m;  100 


Max  z  (m) ;  Nz/meter  =  32 


Figure  15.  Plots  of  Varying  Maximum  z 
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#  2  poims/meter,  iMax  =  1 


Figure  16.  Plots  of  Ez  and  for  Varying  Az 


Two  special  cases  exist,  both  of  which  are  caused  by  rapid  fluctuations  and  large 


estimates  obtained  when  R  becomes  small,  and 


R" 


becomes  large.  The  first  case  occurs 


throughout  the  region  of  interest  when  Ap=p2-pi  becomes  small.  This  causes  the 
calculated  fields  to  be  very  large  and  locally  driven,  even  when  these  calculated  fields  are 
well  behaved  for  a  P2  only  a  few  more  Ap  lengths  away.  Decreasing  the  size  of  Az  is 
required  to  compensate  for  this  effect  as  seen  in  Figure  17.  The  second  case  occurs  when 
both  Pi  and  P2  are  small  in  general.  This  means  that  fields  close  to  the  source  are  being 
estimated.  These  fields  have  rapid  variation  in  the  z  direction  and  have  large  magnitudes 
as  well.  Additionally,  since  pi  is  small  and  pz  is  small,  pi-p2  will  also  be  small,  thus 
exacerbating  the  error  by  including  large  calculated  fields.  Figure  18  shows  that 
decreasing  the  size  of  Az  again  helps  to  counteract  this  effect. 


3.  Varying  the  Number  of  Thin  Ring  Segments 


Varying  the  number  of  segments  (N^)  on  the  ring  leads  to  changes  in  the  segment 
size.  The  segment  size  is  equal  to  PiA(})i. 

PM  =  ^  (82) 

As  increases,  the  length  of  the  segment  size  decreases  and  a  more  accurate  calculation 
for  Un,  Vn,  Wn,  and  the  fields  is  produced.  As  an  initial  estimate,  the  segment  size  should 
be  small  compared  to  of  th®  testing  frequency.  A  more  conservative  limit  on  segment 
size  would  be  to  make  it  less  than  or  equal  to  the  size  of  Az.  Alternatively,  the  segment 
size  can  be  limited  to  between  a  fifth  and  a  tenth  of  the  size  of  Ap=p2-pi.  This  limit 
provides  insight  into  why  segments  should  be  at  least  this  small  if  not  smaller.  As  p2  gets 
closer  to  pi,  the  points  most  influenced  by  the  thin  ring  and  by  the  ring  segments  become 
close  to  the  center  of  the  actual  segment  as  shown  in  Figure  19. 
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Log  Least  Squares  Error  (%);  Msfl(H  )  ^  Least  Squares  Error  (%);  Mag(E^) 


Ljog  Least  Squares  Error  (%)  for  Mag(Ep  for  Varying  Nz  Density  and  p2  ;  p1  =  0.5  m;  2C» 


p2(m) 

#  z  points/meter;  zMax  =  4 


Figure  17.  Case  I;  Mid  range  pi,  with  Small  Ap  and  Moderate  Ap 
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Log  Ua..  Square,  ^  ^  Log  Uas.  Squares  E™r  (%);  Mag(E.) 


Log  Least  Squares  Error  (%)  tor  Mag(Ej)  for  Varying  Nz  Derrsity  and  p2 ;  pi  =  0.01  m;  N^=  100 


#  z  points/meter;  zMax  =  0.5 


Figure  18.  Case  H;  pi  Close  to  Zero,  with  Small  Ap  and  Moderate  Ap 


Large  Segment  Area 


Large  Segment  Area 

Figure  19.  Effect  of  Segment  Size  on  Error 

If  the  overall  size  is  not  reduced,  then  that  segment’s  calculated  value  is  not  a 
good  estimate  for  the  field  contributions  it  represents.  Figure  20  illustrates  the  effects  of 
small  versus  large  N^,  and  large  segment  size  versus  small  segment  size,  respectively. 
Clearly,  for  small  N^,,  the  calculated  fields  are  a  poor  estimate.  Increasing  N,],  quickly 
improves  estimates.  For  further  calculations  the  segment  length  will  be  ~Az  for  small 
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Az.  This  gives  N,  =  for  small  Az.  For  all  calculations  used  in  the  plots,  Na  is  set 

Az 

to  — —  or  100,  whichever  is  the  larger  of  the  two. 
tsz 

4.  Overall  Trends 

From  a  review  of  the  previous  figures,  several  other  interesting  features  can  be 
discussed.  First,  the  calculated  fields  all  have  errors  close  to  the  z  truncation  points. 
This  is  clearly  shown  on  the  plots  of  Relative  Error.  This  is  caused  by  not  including  the 
influence  on  pi  at  points  |  z  I  >  z  cutoff.  By  truncating  fields  greater  than  a  certain  z 
maximum,  error  is  introduced.  The  error  at  the  edges  of  the  calculated  field  is  due  to  the 
proximity  of  field  points  on  pi  that  are  not  included  when  making  an  estimation  of  pa 
fields. 

Another  less  obvious  effect  occurs  when  increasing  p2  from  a  fixed  pi.  After  an 
initial  dip  in  error,  LSE  plots  show  error  slowly  increasing  again  as  P2  increases. 
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Leas!  Squares  Error  (%);  Mag(E^) 
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IV.  BACK  PROPAGATION  OF  FIELDS 


In  forward  propagation  of  fields,  the  contributions  of  thin  rings  of  electric  and 
magnetic  equivalent  currents  determined  from  the  fields  at  the  pi  cylindrical  surface  were 
summed  to  produce  an  estimate  of  the  exact  fields  at  the  p2  cylindrical  surface.  This 
process  was  described  mathematically  in  the  form  of  vector-matrix  multiplication. 


'■EH 


'HH . 


= M4 


'PzCALCUIATED 


(83) 


Algebraically,  to  obtain  the  fields  Ez  and  Hj,  at  pi  from  knowledge  of  the  calculated 
fields  at  p2,  it  appears  that  multiplying  both  sides  of  the  equation  by  the  inverse  of  the 
transfer  matrix  T  will  yield  the  correct  answer. 


(84) 

(85) 


Figure  21  shows  the  effect  of  T*  on  the  calculated  Ez  and  H4,  fields  of  p2.  The  reformed 
Ez  and  Hj,  fields  of  pi  tend  to  be  very  close  to  the  exact  ones.  In  fact,  the  difference  is 
within  machine  accuracy  in  some  areas.  This  is  expected  since  the  T  matrix  actually 
formed  the  calculated  fields. 
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^  ^^.5 Relative  Error  in  Exact  to  Num  Soln  A  for  Mag  INV  E^  {p  1)  from  p  2  Calc  Total  Reids  and  InvTT 


z(p1)  (m);  Density  =  24 ;  p  1  =  0.2  (m);  p  2  =  0.3  (m);  =  100 


^  Relative  Error  in  Exact  to  Num  Soln  A  for  Mag  INV  1)  from  p  2  Calc  Total  Fields  and  frivTT 


2(p1)  (m):  Density  =  24 ;  p  1  =  0.2  (m);  p  2  =  0.3  (m);  *  100 


Figure  21.  Effect  of  T'^  on  the  Calculated  Ez  and  Fields  of  P2. 
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A.  EFFECT  OF  APPLYING  INVERSE  T  MATRIX  TO  THE  EXACT  FIELDS 

ATp2 

When  the  calculated  is  replaced  by  the  exact  ,  a  well  behaved 

T*  would  yield  newly  calculated  inverse  fields  close  to  the  source  fields  at  pi.  If  the 
difference  between  the  calculated  and  exact  fields  at  p2  is  small,  then  applying  a  well 
behaved  T'^  would  yield  a  newly  calculated  field  whose  values  would  be  close  to  the 
exact  pi  field  values.  Mathematically,  the  following  relationship  holds  true  for  a  well 
behaved  and  Sjnvt  is  small. 

(86) 

Unfortunately,  T'^  is  not  well  behaved  in  this  case.  Inverse  T  is  extremely  sensitive  to 
small  changes  in  the  fields  at  P2  used  as  its  input.  As  can  be  seen  from  Figure  22,  even  a 
small  change  of  5a  between  the  calculated  and  exact  fields  at  P2  will  cause  large 
differences  between  the  exact  and  inverted  fields  at  pi.  At  first  glance,  it  seems  that  the 
best  way  to  get  better  results  is  to  make  the  small  change  in  6a  even  smaller.  Figures  23 
and  24  show  the  pointwise  relative  error  of  the  magnitude  of  the  inverted  fields  compared 
to  the  exact  fields  at  pi.  Two  cases  are  shown  in  each  figure,  one  for  a  segment  density 
of  eight  per  meter  and  the  other  for  a  segment  density  of  sixteen  per  meter.  The  LSE  for 
the  inverted  pi  fields  and  the  calculated  p2  fields  is  also  shown  to  give  a  measure  for  the 
overall  improvement  increasing  segment  density  provides.  Also  the  maximum  5a 
between  calcualted  and  exact  fields  at  p2  is  provided. 

Part  of  the  reason  for  the  poor  improvement  in  field  esstimation  lies  with  T  and 
T'\  They  are  classically  ill-conditioned  matrices.  Small  changes  at  the  input  of  an  ill- 
conditioned  matrix  cause  large  changes  at  the  output.  One  measure  of  this  effect  for  a 
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matrix  is  called  the  condition  number.  When  condition  numbers  are  very  large,  matrices 
are  considered  ill-conditioned.  Additionally,  the  transfer  matrix  becomes  physically 
larger  when  either  more  terms  are  used  to  determine  p2  fields  or  when  more  p2  field 
points  are  desired.  Therefore,  as  the  accuracy  of  the  calculated  fields  improves,  more 
points  have  been  used  and  T  becomes  larger  and  larger.  If  T  is  an  ill-conditioned  matrix. 
Figure  25  shows  that  larger  versions  of  T  are  even  more  ill-conditioned.  The  growing 
condition  number  is  one  reason  that  decreasing  5^  has  only  a  small  improvement  in  the 
accuracy  of  the  calculated  inverse  fields.  Therefore,  in  order  to  get  a  better,  smaller  5a  , 
the  size  of  T  must  be  increased,  but  in  increasing  T,  T*  is  increased  and  becomes  even 
more  ill-conditioned,  countering  the  effects  of  a  smaller  5a. 
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(p  1 )  and  INV  (p  1 )  H^(p  1 )  and  INV  H^(p  1 ) 


-2  -1.5  -1  -0.5  0  0.5  1  15  2 

z(p1)  (m);  Density  =  16;  p  1  =  0.2  (m);  p  2  =  0.3  (m);  =  100 


-2  -15  -1  -0.5  0  0.5  1  15  2 

z(p1)  (m);  Density  N  =  16;  p  1  =  0.2  (m);  p  2  =  0.3  (m);  N.  =  100 
z  ^ 


Figure  22.  T*  Sensitivity. 
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Relative  Error  in  Exact  to  Num  Soln  A  for  Mag  INV  E^  (p  1)  from  p  2  Exact  Fields  and  InvTT 
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Relative  Error  in  Exact  to  Num  Soln  A  for  Mag  INV  E^  (p  1)  from  p  2  Exact  Fields  and  InvTT 
-  3r 

N 

UJ 


Leas 

t  Squares 

1 

Mag  INV 

£2  =  21c 

-  n  -igoTC 

.558 

5.  fo 

A 

Mag  E^( 

p2)=00 

m 

m 

■ 

1 

■ 

■ 

■ 

■ 

■ 

A  A 

A  A  A 

[7 

I 

A  A  A 

Mfl 

vA  Aa 

^yv 

lUl 

^  V 

s/ V  V 

TY^ 

vyv^ 

-2  -1.5  -1  -0.5  0  0.5  1  1.5  2 
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Figure  23.  Effect  of  Decreasing  5a  on  Error  in  Electric  Fields 
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-2  -1.5  -1  -0.5  0  0.5  1  1.5  2 


z(p1)  (m):  Density  =  8;  p  1  =  0.2  (m);  p  2  =  0.3  (m);  =  100 


z(p1)  (m);  Density  =  16;  p  1  =  0.2  (m);  p  2  =  0.3  (m);  =  100 


Figure  24.  Effect  of  Decreasing  5^  on  Error  in  Magnetic  Fields 
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Log  Condition  Number  of  INV  TT  for  Varying  Nz  Density 


Figure  25.  Condition  Numbers  of  the  Transfer  Matrix 
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B.  USING  AN  OVERDETERMINED  INVERSE  T  MATRIX 


In  an  attempt  to  improve  the  accuracy  of  the  inverted  fields,  the  inverse  T  matrix 
is  constructed  to  be  over  determined.  The  fields  at  P2  were  determined  at  more  points 
than  the  fields  at  pi.  This  leads  to  a  “long  and  skirmy”  T  matrix  as  indicated  by  the 
following  equation. 

N>M 

Since  T  is  no  longer  a  square  matrix,  the  inverse  of  T  is  not  quite  as  straight  forward. 
Instead  of  the  standard  A'^  of  a  nonsingular  square  matrix,  the  pseudoinverse,  A"^,  is 
calculated  using  the  SVD  technique  shown  in  equation  (88)  where  U  and  V  are  unitary 
matrices  and  S  is  a  diagonal  matrix.  The  pseudoinverse  is  calculated  in  equation  89  from 
the  inverse  of  S,  S^. 

0 

A=USV'  S= 

0 

Jl_ 

•^1 

Sj 

0 

Note,  when  A  is  2Mx2N,  with  N  >  M,  then  A"^  is  2Nx2M  and  AA‘^=I(2Mx2n^-  By  using 
the  SVD  to  find  the  pseudoinverse,  the  matrix  T"^  will  find  the  least  squares  estimate  of 
the  Pi  Inverse  Fields  from  the  over  determined  inputs  [Ref  9].  Figure  26  shows  the  LSE 
inverted  fields  at  pi  when  increasing  the  number  of  P2  segments  by  multiples  of  the 


A^=VS‘tf'  S*  = 
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Least  Squares  Error  (%);  M|ig(INV  H  Least  Squares 


number  of  pi  segments.  After  the  ratio  of  p2  to  pi  segment  densities  increases  beyond 
three,  very  little  improvement  in  the  inverted  fields’  LSE  is  observed.  In  fact,  the  LSE 
degrades  as  the  ratio  increases  from  one  to  three. 


C.  IMPROVING  INTEGRATION  ESTIMATE  ON  pi  WHILE  MAINTAINING 


A  SMALLER  T  MATRIX 


The  next  attempt  at  improving  the  calculated  inverse  field  is  to  maintain  a  smaller 
T  matrix  size  while  attempting  to  get  some  benefit  from  more  integration  accuracy.  In 
order  to  do  this,  a  new  matrix  D  is  defined  such  that  it  is  a  2CN  x  2M  matrix  created 
from  a  z  distribution  as  shown  in  Figure  27.  Then  bands  containing  an  odd  number,  C,  of 
the  matrix’s  rows  are  added  together  and  divided  by  the  number  of  rows  in  the  band,  C. 
This  new  averaged  integration  row  is  then  used  as  a  row  in  T,  a  2N  x  2M  matrix,  where  it 
is  multiplied  by  the  field  point  corresponding  to  the  original  D  matrix  band’s  middle  row. 


r 


C  points  of 
the  D  transfer 
matrix 


point  used  for  Field  Evaluation 
^^^^ntegration  Contributions  Averaged  Together 


(3) 


z  point  used  for  Field  Evaluation 
^^"^ntegration  Contributions  Averaged  Together 


Figure  27.  New  Approximate  Integration  Row  Used  in  T. 
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Thus,  \p\cf,^2M  becomes  by  averaging  bands  of  C  rows,  where  C  is  and  odd 

integer.  Figure  28  shows  the  LSE  of  inverted  fields  at  pi  for  an  increasing  ratio  of  pi 
segment  density  used  to  improve  the  integration  accuracy  to  the  pi  segment  density.  As 
the  ratio  increases  beyond  three,  the  inverted  field  LSE  changes  only  slightly. 


D.  A  CLOSER  LOOK  AT  REAL  AND  IMAGINARY  FIELD  COMPONENTS 


Figures  29  and  30  show  the  Real  and  Imaginary  E-field  components  for  the  exact 
Pi  E-fields  and  the  pi  E-fields  calculated  using  the  matrix  T  inversion  process  for  several 
Pn  P2  pairs.  Figures  31  and  32  show  the  Real  and  hnaginaiy  H-field  components  for  pi 
H-fields.  It  is  striking  that  in  all  cases,  the  Real  Inverted  E-fields  and  the  Imaginary 
Inverted  H-fields  are  essentially  zero. 
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Least  Squares  Error  (%)  for  Mag{lNV  E^)  for  Varying  integration  Nz  Density 


1  3  5  7  9  11  13  15 

Integration  Density  CONST ;  zMax  =  2  m;  pi  =  0.2  m;  p2  =  0.3  m;  N^=  600 


Least  Squares  Error  {%)  for  Mag(jNV  H^)  for  Varying  integration  Nz  Density 


integration  Density  CONST;  zMax  =  2  m;  pi  =  0.2  m;  p2  =  0.3  m;  N^=  600 


Figure  28.  Flat  Weighted  Integration — ^Inverse  Fields 


Z(p1)  (m);  Density  N^,  =  32;  p  1  =  0.2  (m);  p  2  =  0.3  (m);  =  100 


Figure  29.  Back-Propagated  and  Exact  Electric  Field  Components,  Case  1 


Exact  E^(p1)  and  Num  Soln  lor  INV  E^(p1) 


2(p1)(in);  Density  Njj  =  32;  p  1  =  0.5  (m);  p2  =  0.6{m);  N^=  110 


Figure  30.  Back-Propagated  and  Exact  Electric  Field  Components,  Case  2 
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-2  *1.5  -1  -0.5  0  0.5  1  1.5  2 

z(p1)  (m);  Density  N_  =  32;  p  1  =  0.2  (m);  p  2  =  0.3  (m);  N.  =  100 

z  9 

Figure  31.  Back-Progagated  and  Exact  Magnetic  Field  Components,  Case  1 


-4  .3  -2  -1  0  1  2  3  4 

z(p1)  (m);  Density  N  =  32;  p  1  =  0.5  (m);  p2  =  0.6  (m);  N,.=  110 

Figure  32.  Back-Propagaged  and  Exact  Magnetic  Field  Components,  Case  2 
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The  inverted  fields  seem  to  be  composed  of  only  the  imaginary  E-field  and  the 
Real  H-field.  To  investigate  the  cause  of  this,  the  inverted  fields  at  pi  were  propagated 
back  to  p2  to  validate  the  T  matrix.  As  expected  from  a  transform,  inverse  transform  pair 
of  operations,  the  inverted  fields  at  pi  recreated  the  calculated  fields  at  P2.  Thus,  using 
the  T  matrix,  there  are  at  least  two  field  distributions  on  pi  that  create  the  fields  at  P2. 
Upon  further  investigation,  various  combinations  of  the  Real  and  Imaginary  components 
of  the  exact  fields  on  the  pi  cylindrical  surface  were  propagated  out  to  the  P2  cylindrical 
surface  using  the  T  matrix.  Figures  33  through  36  show  that  the  combination  of  the 
Imaginary  E-field  and  Real  H-field  components  of  the  pi  field  distribution  recreated  the 
desired  fields  at  P2.  The  two  field  components  that  formed  the  calculated  P2  field 
distribution  are  also  the  only  non-zero  fields  produced  when  T"’  is  applied  to  the  exact  pi 
fields.  Figures  35  and  36  show  that  the  combination  of  Real  E-field  and  Imaginary  H- 
field  components  of  the  pi  field  distributions  seemed  to  produce  zero  fields  at  pz-  This 
led  to  the  discovery  that  the  real  pi  E-field  and  the  imaginary  pi  H-field  are 
complementary  solutions  to  Maxwell’s  equations  and  are  in  the  null  space  of  the 
transform  matrix  T  operator,  pief  10]  Conversely  this  implies  that  there  is  a  particular 
solution  to  Maxwell’s  equations  that  produces  the  fields  at  pz.  [Ref  10] 
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2(p2)  (m);  Density  =  32;  p  1  =  0.2  (m);  p  2  =  0.3  (m);  =  100 

Figure  33.  Exact  and  Particular  Solution  Component  Calculated  Magnetic  Fields 


Exact  E^(p2)  and  Particular  Soln  for 


2(p2)  (m);  Density  =  32;  p  1  =  0.2  (m);  p  2  =  0.3  (m);  =  100 


Figure  34.  Exact  and  Particular  Solution  Component  Calculated  Electric  Fields 


Exact  E^(p2)  and  Null  Soln  for  E^(p2) 


z{p2)  (m):  Density  N^  =  32;  p  1  =  0.2  (m);  p  2  =  0.3  (m);  =  100 


Figure  35.  Exact  and  Null  Solution  Component  Calculated  Electric  Fields 


Exact  H  (p2)  and  Null  Soln  for  H  (p2) 
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Figure  36.  Exact  and  Null  Solution  Component  Calculated  Magnetic  Fields 
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In  order  to  better  understand  this.  Maxwell’s  equations  are  examined  for  a  strictly 
real  current  and  for  E  and  H-fields  broken  into  real  and  imaginary  components.  While 
J  =  J^+  jjg  in  general,  the  assumed  dipole  current  distribution  is  strictly  real,  J  =  J^. 


Ea  =  Ea2  +  jEai 

(90) 

Ha  =  Hai  +  jH  a2 

(91) 

The  real  and  imaginary  parts  of  Ea  and  Ha  are  indexed  such  that  paired  component 

solutions  have  identical  indices.  Substituting  equations  (90)  and  (91)  into  Maxwell’s 

equations  illustrates  the  respective  field  components  relations. 

V  X  ^a2  +  y£'^i)=  -jco/^Ai  + 

(92) 

Vx(i7^i  +  77/^2)=  ja>£^A2  +  jEAi)+J  A 

(93) 

These  two  equations  are  then  split  into  four  by  grouping  real  and  imaginary  terms. 

V  X  {jEai)=  -JCOJU^Al) 

(94) 

Vx^^i]=  jg}s{jEa\)+Ja 

(95) 

V  X  (£'^2)=  ~jG)fl{jHA2) 

(96) 

Vx  (777.42)=  ja)6^A2)-^0 

(97) 

Thus  Eai  and  Hai  in  equations  (94)  and  (95)  solve  the  particular  solution  for  Maxwell’s 
equations  with  a  real  source  current  J a  .  Similarly,  Eai  and  Ha2  in  equations  (96)  and 
(97)  solve  the  complimentary  or  source  free  solution  for  Maxwell’s  equations.  [Ref  10] 

The  fields  partitioned  into  particular  and  complimentary  solutions  of  Maxwell’s 
equations  can  then  be  transformed  into  equivalent  electric  and  magnetic  current 
distributions  by  the  surface  equivalence  theorem.  The  equivelent  surface  currents  are 
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divided  into  those  created  by  the  particular  fields  and  those  created  by  the  complimentary 
fields  as  shown  in  Figure  37. 


Figure  37.  Equivalent  Siuface  Currents  for  Partioned  Field  Solution 
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From  the  Surface  Equivalence  Theorem,  one  can  now  find  the  influence  of  just  the 
particular  solution  fields  and  propagated  fields  in  Region  2  outside  the  equivalent 
surface.  A  similar  analysis  can  be  done  for  the  complimentary  solution  fields. 


J  SA  and  Msa 


generate 


E  =  0 

E  =  Ea2  +  JEai 
H=HAX+jHA2 


Regionl 

Region2 


JsAi  = 

Msai  =  JEax  X  n 


generate 


E  =  Ea2 
H  =  jHA2 
E  =  Ea2  +  JEaX 
H  =  Max  +  jHA2 


Regionl 

Region2 


J  SA2  —  fix  Ha2 
MsA2  -  j  Ea2  X  h 


generate 


E  =  ~Ea2 
H  =  -JHa2 
E  =  0 
H  =  0 


Regionl 

Region2 


[Ref  10] 


These  results  match  the  results  from  using  the  outward  propagation  transform 
matrix  T  and  the  particular  and  complimentary  fields.  The  above  equations  suggest  that 
using  the  transform  matrix  T  on  particular  solution  fields  will  give  the  complementary 
fields  inside  the  surface.  Put  another  way,  use  the  T  matrix  on  fields  generated  at  P2  to 


get  some  informaiton  about  fields  at  pi.  Figures  38  through  41  confirm  this  result. 
Figures  38  and  39  show  the  contributions  to  the  fields  inside  pi  from  integrating  the 
equivalent  currents  generated  at  pi  by  the  particular  solutions  to  Maxwell’s  equations. 
Inside  pi,  the  particular  solution  at  pi  (the  Real  H-field  and  the  Imaginary  E-field)  create 
the  null  solutions  inside  pi  (the  Imaginary  H-field  and  the  Real  E-field).  Figures  40  and 
41  show  the  contributions  to  the  fields  inside  pi  from  integrating  the  equivalent  currents 
generated  at  pi  by  the  null  solutions  to  Maxwell’s  equations.  Inside  pi,  the  null  solutions 
at  Pi  create  the  inverse  of  the  null  solutions  inside  pi.  Thus,  the  total  field  from  null  and 
particular  solutions  will  cancel  inside  pi.  Unfortunately,  this  method  of  back  propagating 
a  portion  of  the  field  gives  only  complimentary  field  information.  Therefore,  the  inverse 
T  method  is  still  needed  for  back  propagation  of  the  particular  solution. 
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H  (p  2)  and  Num  Soln  H.(p  2) 


Exact  Particular  Soln  for  E  (f^) 
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Figure  38.  Exact  and  Integrated  Particular  Solution  Electric  Fields 


z(p2)  (m);  Density  N  =  32;  p  1  =  0.3  (m);  p  2  =  0.2  (m);  N,  =  100 

Z  7 

Figure  39.  Exact  and  Integrated  Particular  Solution  Magnetic  Fields 
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Exact  Ej,(p2)  and  Null  Soln  for  E^(p2) 


z(p2)  (m);  Density  =  32;  p  1  =  0.3  (m);  p  2  =  0.2  (m);  =  100 


Figure  40.  Exact  and  Integrated  Null  Solution  Electric  Fields 


z(p2)  (m):  Density  N^  =  32;  p  1  =  0.3  (m^  p  2  =  0.2  (m);  =  100 

Figure  41.  Exact  and  Integrated  Null  Solution  Magnetic  Fields 
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V.  CONCLUSION 


As  the  distance  from  a  source  increases,  the  ability  to  resolve  rapid  field 
variations  close  to  the  source  degrades.  Once  in  the  far  field,  two  sources  can  only  be 
distinguished  if  they  are  greater  than  one  half  wavelength  apart.  If  measurements  are 
made  on  the  surface  of  the  object,  the  physical  presence  of  the  probe  disrupts  the  fields  it 
is  attempting  to  sample.  Moving  a  short  distance  away  from  the  source  into  the  near  field 
provides  a  compromise  between  probe  interference  and  resolution  capability.  The  near 
field  probe  will  no  longer  heavily  disturb  the  fields  and  will  be  able  to  measure  fields 
capable  of  providing  better  than  one  half  wavelength  localization  of  somces. 

The  techniques  used  to  analyze  back  propagation  do  not  give  a  perfect 
reconstruction  of  fields  closer  to  the  source.  In  fact,  some  components  of  the  calculated 
inverse  fields  are  essentially  zero  and  provide  no  useful  information  about  the  location  of 
the  source  whatsoever. 

The  reason  that  the  inversion  technique  as  used  in  this  thesis  does  not  give  better 
answers  is  threefold.  At  the  level  of  analyzing  the  transfer  matrix  structure,  the  matrices 
used  are  very  ill-conditioned.  In  fact,  as  the  matrices  become  more  accurate — ^more  z 
range  with  finer  integrations,  they  should  become  “perfectly”  ill-conditioned.  They 
transform  non-zero  vectors  to  zero  exactly.  The  physical  nature  of  the  fields  innately 
creates  ill-conditioned  matrices.  As  more  points  and  finer  resolution  are  used  to  decrease 
error  in  the  calculated  fields  at  pz,  the  transfer  matrix  thus  produced  becomes  more  and 
more  ill-conditioned.  Inverting  the  ill-conditioned  matrix  creates  another  ill-conditioned 
transfer  matrix  for  determining  the  source  fields  at  pi.  The  matrix  is  so  ill-conditioned 
that  even  a  very  small  difference  between  the  calculated  and  exact  field  at  pz  causes  large 
error  in  the  calculated  inverse  field. 

At  an  intermediate  level,  an  actual  radiator  with  a  cylindrical  surface  at  pi  may 
produce  evanescent  fields  as  well  as  the  fields  identical  to  those  from  a  real  source 
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current  located  inside  pi.  These  evanescent  fields  would  have  a  small  but  non-zero 
contribution  to  the  fields  measured  at  p2.  For  equivalent  currents  at  pi,  no  actual 
evanescent  fields  exist.  The  equivalent  currents  are  used  to  make  numerical  estimates  of 
the  fields  at  p2 — ^the  calculated  P2  fields.  When  inverting  the  exact  fields  at  P2  to  estimate 
the  fields  at  pi,  the  small  differences  between  the  calculated  p2  fields  and  the  exact  pz 
fields  could  be  created  by  an  evanescent  field  at  pi.  Thus,  this  technique  has  the 
potential  to  produce  non-physical  evanescent  fields  as  well  as  the  actual  fields  on  pi. 

At  the  highest  level,  there  are  the  cases  of  non-radiating  fields  [Ref  8  and  10]. 
The  fields  have  a  mathematical  structure  such  that  they  can  be  large  at  pi  and  exactly 
zero  at  pz  and  at  all  points  greater  than  pj.  These  fields  are  the  homogeneous  solution  to 
the  forward  transfer  matrix  and  can  be  considered  to  be  in  the  null  space  of  the  forward 
transfer  matrix.  These  fields  can  exist  mathematically,  but  discretizing  the  problem  may 
cause  the  inverse  transfer  matrix  to  reconstruct  the  actual  along  with  some  part  of  a  non¬ 
radiating  or  evanescent  field. 

Therefore,  by  its  very  nature  the  inverse  source  problem  as  handled  in  this 
technique  is  not  only  ill-conditioned,  but  is  perfectly  ill-conditioned,  since  null  or  very 
small  fields  at  pz  can  exist  due  to  large  non-physical  fields  at  p^.  There  will  then  be  the 
actual  inverse  field  plus  any  number  of  non-radiating  and  evanescent  fields  that  can  exist 
at  Pi  and  still  give  the  measured  exact  fields  at  pz- 

The  technique  was  chosen  for  its  ability  to  be  expanded  to  any  axisymmetric 
problem,  not  just  cylindrical  symmetries,  in  a  maimer  that  did  not  exclusively  depend  on 
a  separable  geometry  and  specialized  Hankel  functions. 

There  are  several  areas  for  further  study  using  this  technique.  The  first  is  a  study 
of  modifications  to  the  transfer  function  in  order  to  eliminate  its  null  space  in  the 
inversion  process.  The  second  is  examining  methods  for  decorrelating  the  terms  in  the  T 
matrix  to  simplify  inversion.  The  third  is  wavelet  analysis  of  improved  transfer  matrices, 
specifically  focusing  on  a  wavelet’s  joint  space/spatial  scale  localization  abilities.  The 
fourth  is  applying  this  technique  to  other  axisymmetric  fields  and  sources  (such  as  a 
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spherical  radiator,  conical  radiator,  and  other  sources  found  on  a  surface  of  revolution). 
Finally,  this  technique  can  be  used  to  find  the  actual  surface  currents  on  metal  radiators 
and  scatterers. 
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APPENDIX  A.  METHOD  B-MAXWELL’S  CURL  EQUATIONS  APPLIED 
TWICE  TO  ORIGINAL  VECTOR  POTENTIALS  AND  METHOD  C-SOLVING 

ae7  and  ah  directly  from  vector  potentials  by  an 

ALTERNATIVE  METHOD 


A.  Method  B:  Maxwell’s  Curl  Equations  Applied  Twice  to  Original  Vector 
Potentials 

In  the  previous  section,  AH^  and  Af"^  were  determined  from  Aff^  and  AH^ 
respectively.  In  this  method,  AH^  and  AE^  are  found  directly  without  the  intermediate 
step  of  completely  determining  AE^  and  AH^ .  Since 


E^  =-^Vx/7  = 


1 


then  E. 


Vx 


J(OS, 

1 


Vx 


■VxA 


A 


Vx 


4^  , 


, ,  - jkR 

Ih.—ds' 
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(100) 


(101) 


The  curl  operations  can  both  be  brought  inside  the  integral  to  give  the  following. 
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Breaking  the  surface  of  the  integration  into  many  thin  rings  gives  the  following. 
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In  section  n.B.2  “The  Equivalent  Electric  Current  Case  and  the  Vector  A  Potential,”  it 
was  found  the  following  was  found. 
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*  p,\jc%-^jkV,-^V,\  (107) 


(aeJ,  =  ^^^[2(14 +yic;,)+(fl^ -%iku, -3U,) 

jco  e,  4;r 

+  AP3{6r3+6y«:14-2iV,)J  (108) 
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Note  that  while  (AE^)^  is  the  same  for  Method  A  and  Method  B,  (AE^)^  is  somewhat 
different.  Similarly,  the  following  is  true. 


—  1  v7  1 

- VxE  = - 
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Once  again,  breaking  the  cylindrical  surface  at  pi  into  thin  rings  will  give  AH . 


(109) 
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(111) 


e-"^(^x(i?x(^cos^i-psin<4))) 


-T+  , 

R^  R^ 


X  (i?  X  ((^cos^j  -  psin 


(112) 


(ah,),  =  +  jkV,)*  (pf  +  pI  +  (z,  - zjXeV,  -  ijkV,  - 3V,) 

+  p,pj^u,  +  3  jkU^  -  +3^5+3  jkW^  -  k^W^\  (1 13) 

This  is  the  same  (AH^)^  as  determined  in  Method  A. 
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B.  Method  C:  Solving  AE^  and  AH  Directly  From  Vector  Potentials  by  an 
Alternative  Method. 

In  the  previous  section  AJ?^  and  AH^  were  solved  directly  from  Maxwell’s  curl 

equations  — V  x  and  Hp  =  — — V  y-Ep  .  An  alternative  method  is  to 

jG)S, 

substitute  =  —V  x  A  and  Ep  =  — ^  V  x  F  and  use  Lorentz  condition. 
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Expand  the  second  terms  of  E^  and  Hp . 
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Break  the  cylindrical  surface  into  thin  rings  to  calculate  Hp  and  £'^  and  substitute  Un 
Vji,  and  Wn  functions. 

(AE^)  is  the  same  as  in  the  previous  two  methods. 
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+  —\p,pj3U,+3jkU^ -k^U^  -3^5 -3jkW^^ew^-V^-jkV^  }  (125) 

oae. 


Note  that  while  (AE^)^  is  the  same,  (AE^)^  and  (A//^)^  are  somewhat  different  in  form 
when  compared  to  the  previous  two  methods. 
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APPENDIX  B.  COMPUTER  SOURCE  CODES  DEVELOPED 


%  TestSa.m 

%  By  M.A.  Morgan  13  June  97 

%  Mod-8a  on  30  May  98  substitutes  ifft  for  fft  and  fft  for  ifft 
%  to  correct  kz  polarity. 

%  Propagates  exact  H_phi,  E__z  and  E_rho  fields  from  a  z-axis 
%  array  of  specified  dipoles  between  rhol  and  rho2  using 
%  cylindrical  wave  function  transfer  functions  of  complex  FFT’s. 

%  Fields  are  plotted  on  a  user  specified  2-D  grid  in  rho  and  z. 

%  An  N-element  uniform  array  of  identical  dipole  elements  is  the 
%  source.  A  sequential  phase  shift  on  element  input  currents 
%  provides  beam  steering  of  the  specified  array. 

%  Dipole  element  currents  have  sinusoidal  form  of  I(z)=Im  sin(h~|zl) 
%  which  has  exact  near- field  integration  results  found  in 
%  Jordan  &  Balmain  pp  333-337  or  Elliott  pp  329-332. 

clear  all;  eta=120*pi; 

f=300;%  f=input (’ Enter  frequency  (MHz):  '); 

L=input (’ Enter  dipole  element  length  L  (meters):  ’); 
h=L/2; 

wl=300/f;  k=2*pi/wl;  k2=k*k;  kh=k*h;  Im=l;  ckh=cos(kh); 

N=input ( 'Enter  number  of  dipoles  in  array:  ’); 
if  N  >  1, 

d=input (* Enter  array  spacing  (meters):  ’); 
phi=input ( 'Enter  far-field  pointing  angle  (deg):  '); 

2d=- (N-1) *d/2:d: (N-1) *d/2;  %  Element  positions 

alpha=-k*2d*cos (pi*phi/180) ;  %  Element  phasing 
I=Im*exp ( j*alpha) ;  %  Input  currents 

else, 

2d=0;  I=Im; 

end 

rhol=input ( 'Enter  rhol  (m)  for  field  calculation  grid:  ’); 
rho2=input (' Enter  rho2  (m)  for  field  calculation  grid:  '); 
nr= input (’ Enter  Number  of  rho-Points:  '); 
rho=linspace (rhol,rho2,nr) ;  %  Row  array 

zl=input ( 'Enter  zl  (m)  for  field  calculation  grid:  *); 
z2=input ( 'Enter  z2  (m)  for  field  calculation  grid:  '); 
nz=input (' Enter  Number  of  z-Points  (64,  128,  256,  ...  etc):  ')/ 
dz= (z2-zl) / (nz-1) ;  zee=(zl:dz:z2) ' ;  %  Column  array 

%  Wavenumber  spectra 

dkz=2*pi/ (nz*dz) ;  kz2=nz*dkz/2;  kzl=-kz2+dkz; 
kz=kzl:dkz:kz2;  kzp=0:dkz:kz2;  nkzp=length (kzp) ; 
nsup=find(kzp  <=  k) ;  nsub=find (kzp  >  k) ; 
krsup=sqrt (k2-kzp (nsup) .*kzp (nsup) ) ; 
if  '^isempty  (nsub) , 

krsub=sqrt (kzp (nsub) .*kzp(nsub)-k2) ;  end 
%  Initializing  Arrays 
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FEl==zeros  (nkzp,  1)  ;  FH1=FE1;  Hp=zeros  (nz^nr) ;  Ez=Hp;  Er=Hp; 
[RhO/Zee]=iaeshgrid(rho,zee) ;  Rho2=Rho.*Rho; 

for  n=l:N  %  Array  element  index 

ZnO=Zee“zd(n) ;  Znl=ZnO~h;  Zn2=ZnO+h; 

R0=sqrt(Rho2  +  2n0.^2);  Rl=sqrt (Rho2  +  Znl.^2);  R2=sqrt (Rho2  + 
Zn2.'"2); 

E0=2*exp(~j*k*R0) ;  El-exp {-j*k*Rl) ;  E2-exp(-j*k*R2); 

Hp(:, :)-Hp(:, :)+j*I(n)* (El  +  E2  -  ckh*E0)./{4*pi*Rho); 

Ez(:, :)=Ez(:, :)-j*30*I(n)*(El./Rl  +  E2*/R2  -  ckh*E0,/R0); 

Er(:, ;)=Er(:, : ) + j*30*I (n) * . , . 

(Znl.*El./Rl  +  Zn2.*E2./R2  -  ckh*2n0 • *E0 . /RO) . /Rho; 

end 

svfl-input ( 'Enter  1  to  Save  E  and  H  Field  Files:  •); 
if  -'isempty  (svfl)  ^  if  svfl  —  1, 

hdr=[f  rhol  rho2  nr  zl  z2  nz]; 
save  hdr.dat  hdr  /ascii 
dum-real (Hp) ;  save  hpr.dat  dum  -ascii 
dum=imag(Hp) ;  save  hpi.dat  dum  -ascii 
dum-real (Ez) ;  save  ezr.dat  dum  -ascii 
dum-imag (Ez) ;  save  ezi.dat  dum  -ascii 
dum-real (Er) ;  save  err.dat  dum  -ascii 
dum-imag (Er) ;  save  eri.dat  dum  -ascii 
end;  end 

Fpl-zeros (nz,nr) ;  Ppl-zeros (nz, 1) ;  Pp2-Ppl;  Pp3=Ppl; 

%  FFT  in  z  of  each  fixed  rho  column 
%  using  ifft  in  z  to  correct  polarity  of  k__z  in  FFT 
Hpf=nz*ifft (Hp)  ;  Ezf=nz*ifft (Ez) ;  Erf=nz*if ft (Er) ; 

Pr— pi*Rho.*real (Ez.*conj (Hp) ) ;  %  Time-average  rho-directed  power/dz 

Pz=+pi*Rho,*real (Er.*conj (Hp) ) ;  %  Time-average  z-directed  power/dz 

Prf— Rho.*real (Ezf .*conj (Hpf ) ) /2;  %  Power  Spectral  Densities  (PSD) 
Pzf=+Rho. *real (Erf .*conj (Hpf) ) /2; 


while  1, 


disp ( ' 

Program 

disp ( ’ 

0  — > 

disp ( ' 

1  — > 

disp ( ' 

2  — > 

disp ( ' 

3  — > 

disp ( ’ 

4  — > 

disp ( ’ 

5  — > 

disp ( ’ 

6  — > 

disp ( ' 

7  — > 

disp ( ’ 

8  — > 

disp ( ’ 

9  — > 

disp ( ' 

’) 

Options :  ' ) 

Stop  Program') 

H_phi  3-D  Plot') 

E_2  3-D  Plot') 

E_rho  3-D  Plot’) 

P_rho  3-D  Plot’) 

P_2  3-D  Plot') 

P-vector  Plot') 

Select  Rhol  (Before  Option  8)') 

Fixed  Propagation  from  Rhol  to  Rho2') 
Propagation  Errors  from  Rhol  through  Rho2 ’ 


nplot-input ( 'Select  Numbered  Program  Option:  ’); 
if  isempty (nplot) I  nplot  >  9  |  nplot  <=0,  break; 


end 


if  nplot— 1, 

clf  reset;  surfl (Zee,Rho, abs (Hp) , [0  0]) 
if  N  >  1, 

title  ( [ '  I  H__{phi }  I  for  N= ' ,  int2str  (N) ,  '  Array;  d=  ’ ,  num2str  (d) ,  .  • . 
’  m;  L=’ /niam2str  (L) ,  ’  m;  phi=’,nxam2str  (phi) , 'deg;  f-'/,.. 
num2str(f),'  MHz']) 
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else, 

title  ( [ '  lH__{phi}  I  for  L=' ,  num2str  (L) ,  ’  m  Dipole  at  f=' ,num2str  (f) , 
'  MHz']) 

end 

xlabelCz  (m)  ' ) ;  ylabel('rho  (m)  ' ) ;  zlabel  (' |H_{phi}  |  (A/m)') 

view (200, 30) ;  v=axis;  axis([zl  z2  rhol  rho2  v(5)  v(6)]);  figure{l) 

hep y=input (' Enter  1  for  Hard  Copy:  ' ) ; 

if  ~isempty (hepy) ,  if  hepy  ==  1,  print;  end;  end 

shading  interp;  view (90, -90) ;  colorbar;  figure(l) 

hep y=input ( 'Enter  1  for  Hard  Copy:  '); 

if  ~isempty (hepy) ,  if  hepy  ==  1,  print;  end;  end 

elf  reset;  surf (Zee, Rho, loglO (abs (Hp) ) ) 
if  N  >  1, 

title ([' |H_{phi} 1  for  N=' , int2str (N) , '  Array;  d=' ,num2str (d) ,  . . . 

'  m;  L=' ,num2str (L) , '  m;  phi=' ,num2str (phi) , 'deg;  f=',... 
num2str(f),'  MHz']) 

else, 

title ( [ ' |H_{phi} I  for  L=' , num2str (L) , '  m  Dipole  at  f=' ,num2str (f) , 
'  MHz'] ) 

end 

xlabel ( ' z  (m)  ' ) ;  ylabel ( ' rho  (m)  ' ) 
v=axis;  view(90,-90) ;  shading  interp 

axis([zl  z2  rhol  rho2  v(5)  v(6)]);  eolorbar;  figure (1) 
hepy=input ( 'Enter  1  for  Hard  Copy:  '); 
if  -isempty (hepy) ,  if  hepy  ==  1,  print;  end;  end 

%  Manual  fftshift  for  1-D  fft  array 
Fpl (1 :nkzp-2,  : ) =abs (Hpf (nkzp+1 :nz, : ) ) ; 

Fpl (nkzp-1 :nz,  : ) =abs (Hpf (1 :nkzp, : ) ) ; 
elf  reset;  surfl (kz, rho, Fpl ' , [0  0]) 
if  N  >  1, 

title ([ 'H_(phi}  Speetrum  for  N=' , int2str (N) , '  Array; 
d=' ,num2str (d) , . . . 

'  m;  L=' ,n\am2str  (L)  ,  '  m;  phi=' ,n\im2str  (phi) , 'deg;  f=',... 
num2str(f),'  MHz']) 

else, 

title ([ 'H_{phi}  Speetrum  for  L=' ,num2str (L) , '  m  Dipole  at  f=',... 
num2str(f),'  MHz']) 

end 

xlabel ('k_z  (rad/m)');  ylabel ('rho  (m) ' ) ;  zlabel CIFFT  [H_{phi}]|') 

view (200, 30) ;  v=axis;  axis([kzl  kz2  rhol  rho2  v(5)  v(6)]);  figure (1) 

hep y=input (' Enter  1  for  Hard  Copy:  '); 

if  ~isempty (hepy) ,  if  hepy  ==  1,  print;  end;  end 

shading  interp;  view(90,-90) ;  eolorbar;  figure(l) 

hep y=input ( 'Enter  1  for  Hard  Copy:  '); 

if  "isempty (hepy) ,  if  hepy  ==  1,  print;  end;  end 

elf  reset;  surf (kz, rho, loglO (Fpl ' + . 01) ) 
if  N  >  1, 

title ([ 'H_{phi}  Speetrim  for  N=' , int2str (N) , '  Array; 
d=' ,num2str (d) ,  . . . 

'  m;  L=' ,num2str (L) , '  m;  phi=' , num2str (phi)  ,  ' deg;  f=', _ 

num2str(f) ,  '  MHz'] ) 

else, 

title ([ 'H_{phi}  Speetrum  for  L=' ,num2str (L) , '  m  Dipole  at  f=',... 
num2str(f),'  MHz']) 
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end 

xlabel ( ’ k^z  (rad/m) ' ) ;  ylabel ( ’ rho  (m) ' ) 
v=axis;  view(90, -90) ;  shading  interp 

axis([k2l  kz2  rhol  rho2  v(5)  v(6)]);  colorbar;  figure (1) 
hep y=input  (* Enter  1  for  Hard  Copy:  M; 
if  '-isempty  (hepy) ,  if  hepy  ==  If  print;  end;  end 
end 

if  nplot==2, 

clf  reset;  surfl (Zee,RhO/ abs (Ez) , [0  0]) 
if  N  >  1, 

title  ( [  ’  |E_2  I  for  N=' ,  int2str  (N) ,  '  Array;  d=’ ,num2str  (d) ,  '  m; 

num2str(L),’  m;  phi=’ ,num2str  (phi) , ’deg;  f=’ ,num2str  (f ) ,  *  MHz’]) 

else, 

title  ( [  ’  I  E_z  1  for  L=’ ,num2str  (L) ,  ’  m  Dipole  at  f=’ ,nuia2str  (f ) ,  ’ 
MHz’] ) 
end 

xlabel  (’z  (m)  ’ ) ;  ylabel  (’rho  (m)  ’ ) ;  zlabel  ( ’  |E_z  |  (V/m)  ’ ) 

view (200, 30) ;  v=axis;  axis([zl  z2  rhol  rho2  v(5)  v(6)]);  figure (1) 

hcpy=input (’ Enter  1  for  Hard  Copy:  ’); 

if  -'isempty  (hepy) ,  if  hepy  ==  1,  print;  end;  end 

shading  interp;  view(90, -90) ;  colorbar;  figure(l) 

hep y=input ( ’Enter  1  for  Hard  Copy:  ’); 

if  -'isempty (hepy) ,  if  hepy  ==  1,  print;  end;  end 

clf  reset;  surf (Zee, Rho, loglO (abs (Ez) ) ) 
if  N  >  1, 

title  ( [  ’  |E_2  I  for  N=’ ,  int2str  (N) ,  ’  Array;  d=’ ,nuia2str  (d) ,  .  •  . 

’  m;  L=’ ,num2str  (L) ,  ’  m;  phi=’ ,  num2str  (phi) , ’deg;  f  =’,... 
num2str (f ) , ’  MHz’]) 

else, 

title  ( [  ’  |E_z  I  for  L=’ ,num2str  (L) ,  ’  m  Dipole  at  f=’ ,nxam2str  (f ) , _ 

’  MHz’]) 

end 

xlabel ( ’ z  (m)  ’ ) ;  ylabel ( ’ rho  (m) ’ ) 
v=axis;  view(90,-90) ;  shading  interp 

axis([zl  z2  rhol  rho2  v(5)  v(6)]);  colorbar;  figure (1) 
hep y=input (’ Enter  1  for  Hard  Copy:  ’); 
if  -'isempty  (hepy) ,  if  hepy  ==  1,  print;  end;  end 

%  Manual  fftshift  for  1-D  fft  array 
F^l (1 :nk2p-2, : ) =abs (Ezf (nkzp+1 :nz,  : )  )  ; 

Fpl (nkzp-1 :nz, : ) =abs (Ezf (1 :nkzp, ; ) ) ; 
clf  reset;  surfl (kz, rho, Fpl ’ , [0  0]) 
if  N  >  1, 

title ([’E_2  Spectrum  for  N=’ , int2str (N) , ’  Array;  d=’ ,num2str (d) , . . . 
’  m;  L=’ ,num2str  (L) ,  ’  m;  phi=’ ,num2str  (phi) , ’deg;  f=’,... 
num2str(f),’  MHz’]) 

else, 

title ([’E_2  Spectrum  for  L=’ ,num2str (L) , ’  m  Dipole  at  f=’,... 
num2str(f),’  MHz’]) 

end 

xlabel  ( ’  k__z  (rad/m)  ’ )  ;  ylabel  ( ’  rho  (m)  ’ )  ;  zlabel  ( ’  |  FFT  tE__z]  |  ’ ) 
view(200, 30)  ;  v=axis;  axis([kzl  kz2  rhol  rho2  v(5)  v(6)]);  figured) 
hep y=input (’ Enter  1  for  Hard  Copy:  ’); 
if  '-isempty  (hepy) ,  if  hepy  ==  1,  print;  end;  end 
shading  interp;  view(90,-90) ;  colorbar;  figure(l) 


82 


hcpy=input (' Enter  1  for  Hard  Copy:  ’); 
if  ~isempty (hcpy) /  if  hcpy  ==  1/  print;  end;  end 

clf  reset;  surf (kz, rho, loglO (Fpl '  +  . 01) ) 
if  N  >  1, 

title (['E_z  Spectrum  for  N=' , int2str (N) ,  '  Array;  d=' ,num2str (d) , . 
'  m;  L=' ,num2str (L) , '  m;  phi=' , num2str (phi) , ’ deg;  f=',... 
num2str(f), '  MHz']) 

else, 

title (['E_z  Spectrum  for  L= ' , num2str (L) ,  '  m  Dipole  at  f=',... 
num2str(f) , '  MHz'] ) 

end 

xlabel ( ' k_z  (rad/m) ' ) ;  ylabel ( | rho  (m) ' ) 
v=axis;  view(90,-90) ;  shading  interp 

axis{[kzl  kz2  rhol  rho2  v(5)  v(6)]);  colorbar;  figure (1) 
hcpy=input (' Enter  1  for  Hard  Copy:  '); 
if  ~isempty (hcpy) ,  if  hcpy  ==  1,  print;  end;  end 
end 

if  nplot==3, 

clf  reset;  surf 1 (Zee, Rho, abs (Er) ,  [0  0]) 
if  N  >  1, 

title ( [ ' lE_{rho} I  for  N= ' , int2str (N) , '  Array;  d=' , num2str (d) , . . . 

'  m;  L=' ,num2str (L) , '  m;  phi=' ,num2str (phi) , . . , 

'deg;  f=' ,num2str (f ) , '  MHz']); 

else, 

title  ( [ '  |E_{rho}  I  for  L=' , num2str  (L)  ,  '  m  Dipole  at  f=' ,nxm:i2str  (f) 
'  MHz ' ] ) 

end 

xlabel  Cz  (m)  ' )  ;  ylabel  ('rho  (m)  ' )  ;  zlabel  ( '  |  E_{rho}  |  (V/m)  ' ) 

view(200,30)  ;  v=axis;  axis([zl  z2  rhol  rho2  v(5)  v(6)]);  figured) 

hep y=input ( 'Enter  1  for  Hard  Copy:  '); 

if  ~isempty(hcpy) ,  if  hcpy  ==  1,  print;  end;  end 

shading  interp;  view (90, -90) ;  colorbar;  figure(l) 

hep y=input ( 'Enter  1  for  Hard  Copy:  '); 

if  ~isempty (hcpy) ,  if  hcpy  ==  1,  print;  end;  end 

clf  reset;  surf (Zee, Rho, loglO (abs (Er) ) ) 
if  N  >  1, 

title ( [ ' |E_{rho} I  for  N=' , int2str (N) , '  Array;  d=' ,num2str (d) , . . . 

'  m;  L=' ,nimL2str  (L) ,  '  m;  phi=' ,num2str  (phi) , 'deg;  f=',... 
num2str(f) , '  MHz'] ) 

else, 

title([' |E_{rho} I  for  L=' , num2str (L) ,  '  m  Dipole  at  f=' ,num2str (f ) 
'  MHz' ] ) 

end 

xlabel ( ' z  (m)  ' ) ;  ylabel ( ' rho  (m) ' ) 
v=axis;  view(90,-90) ;  shading  interp 

axis([zl  z2  rhol  rho2  v(5)  v(6)]);  colorbar;  figure (1) 
hcpy=input ( 'Enter  1  for  Hard  Copy:  '); 
if  -isempty (hcpy) ,  if  hcpy  ==  1,  print;  end;  end 

%  Manual  fftshift  for  1-D  fft  array 
Fpl (1 :nkzp-2, : ) =abs (Erf (nkzp+1 :nz,  : ) )  ; 

Fpl (nkzp-1 :nz, : ) =abs (Erf (1 :nkzp, : ) )  ; 
clf  reset;  surfl (kz, rho, Fpl ' , [0  0]) 
if  N  >  1, 
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title ([ ’E_{rho}  Spectrum  for  N=’ , int2str (N) , ’  Array; 
d=' ,num2str (d) , . • • 

’  m;  L=’ ,num2str (L) , ^  m;  phi=' ,num2str (phi) , 'deg; 
num2str (f ) , '  MHz ' ] ) 

else, 

title ([ ’E_{rho}  Spectrum  for  L=' ,num2str (L) , ’  m  Dipole  at  f=',... 
num2str(f),'  MHz']) 

end 

xlabel('k_2  (rad/m)');  ylabel(’rho  (m)  ' ) ;  zlabel('|FFT  [E_{rho}]|’) 

view(200, 30) ;  v^axis;  axis([kzl  kz2  rhol  rho2  v(5)  v(6)]);  figure (1) 

hcpy=input ( 'Enter  1  for  Hard  Copy:  '); 

if  '^isempty  (hcpy) ,  if  hcpy  —  1,  print;  end;  end 

shading  interp;  view {90, -90)  ;  colorbar;  figured) 

hcpy=input ( 'Enter  1  for  Hard  Copy:  '); 

if  --isempty (hcpy) ,  if  hcpy  ==  1,  print;  end;  end 

clf  reset;  surf (kz, rho, loglO (Fpl '  +  .01) ) 
if  N  >  1, 

title ([ 'E_{rho}  Spectrum  for  N=' , int2str (N) , '  Array; 
d=' ,num2str (d) , . • • 

'  m;  L=' ,num2str (L) , '  m;  phi=' ,num2str (phi) , 'deg;  f=’, _ 

num2str(f),'  MHz']) 

else, 

title  ([ 'E__{rho}  Spectrum  for  L=' ,num2str  (L) ,  '  m  Dipole  at  f 
num2str (f ) , ’  MHz ' ] ) 

end 

xlabel  ( '  k__2  (rad/m)  ’ )  ;  ylabel  ( ’  rho  (m)  ’ ) 
v=axis;  view (90, -90) ;  shading  interp 

axis([kzl  k22  rhol  rho2  v(5)  v(6)]);  colorbar;  figure (1) 
hcpy=input ( 'Enter  1  for  Hard  Copy:  '); 
if  --isempty  (hcpy) ,  if  hcpy  ==  1,  print;  end;  end 
end 

if  nplot==4, 

clf  reset;  surf 1 (Zee, Rho, Pr, [0  0]) 
if  N  >  1, 

title ([ 'P_{rho}  for  N=' , int2str (N) , ’  Array;  d=' ,num2str (d) , . . . 

'  m;  L=' ,num2str (L) , '  m;  phi=' ,num2str (phi) , • . . 

'deg;  f=' ,num2str (f ) , '  MHz']); 

else, 

title ([ 'P_{ rho}  for  L=',num2str (L) , '  m  Dipole  at  f=' ,num2str  (f ) , . . 
'  MHz']  ) 

end 

xlabel  ('z  (m)  ' )  ;  ylabel  ('rho  (m)  ' ) ;  zlabel  ( ' P__{rho}  (W/m)  ' ) 

view (200, 30) ;  v=axis;  axis{[zl  z2  rhol  rho2  v(5)  v(6)]);  figure(l) 

hep y=input { 'Enter  1  for  Hard  Copy:  '); 

if  --isempty  (hcpy) ,  if  hcpy  ==  1,  print;  end;  end 

shading  interp;  view (90, -90) ;  colorbar;  figure(l) 

hcpy=input ( 'Enter  1  for  Hard  Copy:  '); 

if  --isempty  (hcpy) ,  if  hcpy  ==  1,  print;  end;  end 

clf  reset;  surf (Zee, Rho, loglO (abs (Pr)  )  ) 
if  N  >  1, 

title ([ 'P_{ rho}  for  N=' , int2str (N) , '  Array;  d=' ,num2str (d) , . .  • 

'  m;  L=',num2str(L) , '  m;  phi=' ,num2str (phi) , 'deg;  f=',..* 
num2str(f),'  MHz']) 

else. 
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title  {[ ’P_{rho}  for  L=’ ,nuia2str  (L) ,  '  m  Dipole  at  f=*  ,niairi2str  (f ) ,  . . . 
’  MHz’]) 

end 

xlabel ( ’ z  (m)  ’ ) /  ylabel ( ’ rho  (m) ’ ) 
v=axis;  view(90,-“90) ;  shading  interp 

axis([zl  z2  rhol  rho2  v(5)  v{6)]);  colorbar;  figure (1) 
hcpy=input ( 'Enter  1  for  Hard  Copy:  ’); 
if  -isempty (hcpy) ,  if  hcpy  ==  1/  print;  end;  end 

%  Manual  fftshift  for  1-D  fft  array 

Fpl (1 :nkzp-2, : ) =Prf (nkzp+1 :nz, : ) ;  Fpl (nkzp-1 rnz, : )=Prf (1 znkzp, : ) ; 
clf  reset;  surfl (kz, rho, Fpl ’ , [0  0]) 
if  N  >  1, 

title ([’P_r  PSD  for  N=’ , int2str (N) , ’  Array;  d=’ ,num2str (d) , , • . 

’  m;  L=' ,nuici2str  (L) ,  ’  m;  phi=’ ,num2str  (phi)  , 'deg;  f=’,... 
nuin2str (f ) ,  ’  MHz’]) 

else, 

title ([’P_r  PSD  for  L=’ ,num2str (L) , '  m  Dipole  at  f=',,.. 
num2str (f ) , ’  MHz’ ] ) 

end 

xlabel ( ’ k_z  (rad/m) ’ ) ;  ylabel { ’rho  (m)  ’ )  ;  zlabel ( ’P_r  PSD’ ) 

view(200, 30)  ;  v=axis;  axis([kzl  k22  rhol  rho2  v(5)  v(6)]);  figured) 

hcpy=input ( ’Enter  1  for  Hard  Copy:  ’); 

if  -isempty  (hcpy) ,  if  hcpy  ===  1,  print;  end;  end 

shading  interp;  view(90, -90) ;  colorbar;  figure(l) 

hcpy=input ( ’Enter  1  for  Hard  Copy:  '); 

if  -isempty  (hcpy) ,  if  hcpy  =-  1,  print;  end;  end 

clf  reset;  surf (kz,rho, loglO (abs (Fpl) ’+.01) ) 
if  N  >  1, 

title ([’P_r  PSD  for  N=’ , int2str (N) , ’  Array;  d=’ ,num2str (d)  ,  .  . . 

’  m;  L=’ ,nuia2str  (L) ,  ’  m;  phi=’ ,num2str  (phi) , ’deg;  f=’,... 
nuia2str  (f ) ,  ’  MHz’]) 

else, 

title  ([’P_r  PSD  for  L=’ ,nuin2str  (L) ,  ’  m  Dipole  at  f=’,... 
num2str(f),’  MHz’]) 

end 

xlabel  ( ’  k_z  (rad/iti)  ’ ) ;  ylabel  ( ’  rho  (la)  ’ ) 
v==axis;  view(90,-90)  ;  shading  interp 

axis([kzl  kz2  rhol  rho2  v(5)  v(6)]);  colorbar;  figure (1) 
hcpy=input ( ’Enter  1  for  Hard  Copy:  ’); 
if  -'isempty  (hcpy) ,  if  hcpy  ==  1,  print;  end;  end 
end 

if  nplot~5, 

clf  reset;  surfl (Zee,Rho,Pz, [0  0]) 
if  N  >  1, 

title([’P_2  for  N=’ ,  int2str  (N) ,  ’  Array;  d==’ ,num2str  (d) ,  • . . 

’  m;  L=’ ,nuiu2str  (L) ,  ’  m;  phi=’ ,num2str  (phi) ,  .  .  . 

’deg;  f=’ ,num2str  (f ) ,  ’  MHz’]); 

else, 

title  ([’P_2  for  L=’ ,nuia2str  (L)  ,  ’  m  Dipole  at  f=’ ,num2str  (f )  ,  .  . . 

’  MHz’]) 

end 

xlabel  ( ’  z  (m)  ' ) ;  ylabel  ( ’  rho  (m)  ’ )  ;  zlabel  ( ’  P_2  (W/m)  ' ) 

view (200, 30) ;  v=axis;  axis([zl  z2  rhol  rho2  v(5)  v{6)]);  figure{l) 

hcpy=input ( ’Enter  1  for  Hard  Copy:  ’); 
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if  '-isempty (hcpy)  /  if  hcpy  ==  1,  print;  end;  end 
shading  interp;  view(90, -90) ;  colorbar;  figure(l) 
hep y=input ( ’Enter  1  for  Hard  Copy:  *); 
if  '-isempty  (hcpy)  /  if  hcpy  ==  print;  end;  end 

clf  reset;  surf (Zee, Rho, loglO (abs (Pz) ) ) 
if  N  >  1, 

title  ([’P_z  for  N=* ,  int2str  (N) ,  ’  Array;  d=’ ,nuia2str  (d) ,  . . . 

’  m;  L=’ ,nuin2str  (L)  ,  '  m;  phi=’ ,nuin2str  (phi) , ’deg;  f=’,... 
nuni2str  (f ) ,  ’  MHz’]) 

else, 

title  ([’P_z  for  L=’ ,nuin2str  (L) ,  ’  m  Dipole  at  f=’ ,n™2str  (f ) ,  . . . 

’  MHz’]) 

end 

xlabel ( ’ z  (m)  ’ ) ;  ylabel ( ’ rho  (m)  ’ ) 
v=axis;  view(90, -90) ;  shading  interp 

axis([zl  z2  rhol  rho2  v(5)  v(6)]);  colorbar;  figure (1) 
hcpy=input ( ’Enter  1  for  Hard  Copy:  ’); 
if  -isempty (hcpy) ,  if  hcpy  ==  1,  print;  end;  end 

%  Manual  fftshift  for  1-D  fft  array 

Fpl  (l:nkzp-2, : )=Pzf (nkzp+l:nz, : ) ;  Fpl (nkzp-l:nz, : ) =Pzf (1 :nkzp, : ) ; 
clf  reset;  surfl (kz, rho, Fpl ’ , [0  0]) 
if  N  >  1, 

title ([’P_z  PSD  for  N=’ , int2str (N) , ’  Array;  d=’ ,num2str (d) , . . . 

’  m;  L=’ ,num2str  (L)  ,  ’  m;  phi=’ ,  num2str  (phi) , ’deg;  f=’, — 
num2str (f ) , ’  MHz’]) 

else, 

title ([’P_z  PSD  for  L=’ ,num2str (L) , ’  m  Dipole  at  f=’, — 
nuia2str  (f ) ,  ’  MHz’]) 

end 

xlabel  (’k_z  (rad/m)’);  ylabel  (’rho  (in)’);  zlabel(’P__z  PSD’) 

view(200,30) ;  v=axis;  axis([kzl  kz2  rhol  rho2  v(5)  v(6)]);  figure (1) 

hep y=input ( ’Enter  1  for  Hard  Copy:  ’); 

if  -isempty (hcpy) ,  if  hcpy  ==  1,  print;  end;  end 

shading  interp;  view (90, -90) ;  colorbar;  figure(l) 

hep y=input (’ Enter  1  for  Hard  Copy:  ’); 

if  -isempty (hcpy) ,  if  hcpy  ==  1,  print;  end;  end 

clf  reset;  surf (kz, rho, loglO (abs (Fpl) . 01) ) 
if  N  >  1, 

title  (['P_z  PSD  for  N=’ ,  int2str  (N) ,  ’  Array;  d=’ ,nuin2str  (d) ,  . . . 

’  m;  L=' ,nuin2str  (L) ,  '  m;  phi=’ ,num2str  (phi) , ’deg;  f=’,.,. 
num2str  (f ) ,  '  MHz’]) 

else, 

title  ([’P_z  PSD  for  L=’ ,nuiti2str  (L) ,  ’  m  Dipole  at  f=', — 
num2str  (f ) ,  ’  MHz’]) 

end 

xlabel ( ’ k_z  (rad/m)  ’ ) ;  ylabel ( ’ rho  (m) ’ ) 
v=axis;  view(90, -90) ;  shading  interp 

axis((kzl  kz2  rhol  rho2  v(5)  v(6)]);  colorbar;  figure (1) 
hep y=input ( ’Enter  1  for  Hard  Copy:  ’); 
if  --isempty  (hcpy) ,  if  hcpy  ~  1,  print;  end;  end 
end 

if  nplot==6, 

npwr=input ( 'Enter  inverse  scaling  factor  for  quiver  plot:  ’) 
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Prmax=max (max (abs (Pr) ) ) ;  Pzmax=max (max (abs (Pz) ) ) ; 

Pmax=min  (Prmax,  Pzmax) ; 

PrO=sign  (Pr)  .  *  (abs  (Pr)  /Pmax)  .  (1/npwr)  ; 

PzO=sign(Pz)  .*  (abs  (Pz)  /Pmax)  (1/npwr)  ; 
clf  reset;  quiver (Rho, Zee, PrO, PzO) 
if  N  >  1, 

title (['P-vector  for  N=’ , int2str (N) , ’  Array;  d= ' , numZstr (d) , . . . 

'  m;  L=',num2str(L) ,  '  m;  phi=' ,  niam2str  (plii) ,  '  deg;  f=',... 
num2str(f)/'  MHz']) 

else, 

title ([ 'P-vector  for  L=' ,num2str (L) , '  m  Dipole  at  f=' , num2str (f ) ,  . . . 
'  MHz'  ] ) 

end 

xlabel ( ' rho  (m) ' ) ;  ylabel ( ' z  (m) ' ) 

axis([rhol  rho2,  zl  z2] )  ;  axis  square;  figured) 

hcpy=input( 'Enter  1  for  Hard  Copy:  d; 

if  ~isempty (hcpy) ,  if  hcpy  ==  1,  print;  end;  end 

end 

if  nplot==7, 

rl=input( 'Enter  Rhol  radius  to  propagate  from:  '); 

%  Selecting  closest  grid  radius 

[rm,  ml]=min(abs (rho-rl) ) ;  rpl=rho(ml); 

%  Forming  denominator  of  radial  transfer  functions 
FEl=zeros (nkzp,  1) ;  FH1=FE1; 

[H,  DH]=dhanlcel  (0,  krsup*rpl)  ;  FEl(nsup)=H;  FHl  (nsup)  =DH; 
if  ~isempty (nsub) , 

[K,  DK]=dmodbes (0, krsub*rpl) ;  FEl(nsub)=K;  FHl (nsub) =DK;  end 

%  Energy  Normalization  for  Optimal  Estimators 
EE1=FE1 ' *FEl/nkzp;  EH1=FH1 ' *FHl/nkzp; 

%  Plotting  Spatial  and  Spectral  Fields  at  rho=rpl 
clf  reset 
subplot (2,1,1), 

plot  (zee, abs (Ez(:,  ml) ) , ' g' ,  zee,  abs (eta*Hp ( :  ,ml) ) ,  '  :r',  . . . 

zee, abs (Er ( : ,ml) ) , ' — c' ) 
if  N  >  1, 

title ( ['Fields  and  Spectra;  Array  N=',int2str(N),'; 
d=' ,num2str (d) ,  . . . 

'm;  L=',num2str (L) , 'm;  phi=' ,num2str (phi) , 'deg;  f=',... 
num2str(f) , 'MHz;  rho=' ,num2str (rpl) , 'm'] ) 

else, 

title ( ['Fields  and  Spectra  for  L=' ,num2str (L) , 'm  Dipole  at  f=',... 
num2str (f ) , 'MHz;  rho=' ,nxmi2str (rpl) , 'm' ] ) 

end 

xlabel ( ' z  (m) ' ) ;  ylabel ( ' V/m' ) 

Ppl(l:nkzp-2)=abs(Ezf (nkzp+l:nz,ml) ) ;  Ppl (nkzp- 
l:nz)=abs (Ezf (1 :nkzp,ml) ) ; 

Pp2  (l:nkzp-2)=eta*abs(Hpf (nkzp+l:nz,ml) ) ;  %  proper  fftshift  operation 
Pp2 (nkzp-1 :nz) =eta*abs (Hpf (1 :nkzp,ml) ) ; 

Pp3(l:nkzp-2)=abs(Erf (nkzp+l:nz,ml) ) ;  Pp3(nkzp- 
1 ;nz) =abs (Erf (1 : nkzp, ml) ) ; 

subplot  (2, 1,2) ,  plot(kz,Ppl,  'g',kz,Pp2,  '  :r',kz,Pp3,  '— c') 
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title  ( ’  |E_2  I  (solid)  ;  eta*  |H__{phi}  |  (short  dash) ;  lE_{rho}  |  (long 
dash) ’ ) 

xlabel ( ’ k_z  (rad/m) ’ ) ;  ylabel ( ’ 1 FFT I ’ ) ;  figure (1) 

^cpy=iuput ( ’Enter  1  for  Hard  Copy:  ’); 
if  -'isempty  (hcpy) ,  if  hcpy  ==  If  print;  end;  end 

clf  reset 

subplot  (2,1,1),  plot  (zee,  Pr  ( :  ,ial) ,  ’  g’ ,  zee,  Pz  ( :  ,inl) ,  *r :  ’ ) 
if  N  >  1, 

title ([ ’Power  &  PSD;  Array  N=’ , int2str (N) ,’ ;  d=’ ,num2str (d) _ _ 

’m;  L=’  ,num2str  (L)  ,  ’la;  phi”’  ,nuin2str  (phi) ,  ’deg;  f=’ ,  , . . 
nuia2str  (f )  ,  ’MHz;  rho=’  ,nuin2str  (rpl) , ’m’  ] ) 

else, 

title  ([’ Fields  and  Spectra  for  L“’ ,nuia2str  (L) , ’m  Dipole  at  f=*,.*, 
nuia2str  (f ) ,  ’MHz;  rho”’  ,nuiri2str  (rpl) , ’m’  ] ) 

end 

xlabel  (’z  (m)  ’ ) ;  ylabel  ( ’W/m’ ) 

%  proper  fftshift  operation 

Ppl  (l:nkzp-2)=Prf  (nk:zp+l:nz,ml) ;  Ppl  (nkzp~l:nz)=Prf  (l:nkzp,ial)  ; 

Pp2 (l:nkzp-2)=Pzf (nkzp+l:nz,ml) ;  Pp2 (nkzp-1 :n2)=Pzf (1 :nkzp,ml)  ; 

subplot (2,1,2),  plot (kz, Ppl, ’g’ , kz, Pp2, ’r : ’ ) 

title (’ Plots :  P_r  (solid)  P_z  (dash)’) 

xlabel  (’k__z  (rad/ia)*);  ylabel  ( ’PSD’ )  ;  figure  (1) 

hcpy=input ( ’Enter  1  for  Hard  Copy:  ’ ) ; 

if  -isempty (hcpy) ,  if  hcpy  ==  1,  print;  end;  end 

end 

if  nplot==8, 

if  isempty(rl),  error (’Enter  Rhol  First’);  end 
r2=input ( ’Enter  Rho2  radius  to  propagate  to:  ’); 

%  Selecting  closest  grid  radius 
[  rm,  ia2  ]  =min  { abs  ( rhO“r2 )  )  ;  rp2=rho  (m2 )  ; 

%  Initialize  Arrays 

TEz=zeros (nz, 1) ;  THp=TEz;  Ppl=TE2;  Pp2=TEz; 

FE2=zeros (nkzp, 1) ;  FH2=FE1; 

%  Forming  numerator  of  radial  transfer  functions 

[H,  DH]=dhankel (0, krsup*rp2) ;  FE2(nsup)=H;  FH2  {nsup)=DH; 

if  -'isempty  (nsub) , 

[K,  DK]=dmodbes  (0,  krsub*rp2)  ;  FE2(nsub)=K;  FH2  (nsub)=DK;  end 
while  2 

CEz=input ( ’Enter  Optimal  Estimator  Constant  for  E_2  (CEz  <  0  to  exit): 

M; 

if  CEz  <  0,  break;  end 

CHp^input ( ’Enter  Optimal  Estimator  Constant  for  H_phi  &  E_rho:  ’); 

%  Computing  Riad’s  Optimal  Deconvolution  Propagation  Transfer  Functions 
TEz  (l:nkzp)=FE2.*conj  (FED  ./ (FEl.*conj  (FED  +  CEz*EED; 
THp(l:nkzp)=FH2,*conj (FHl) ./(FHl.*conj (FHl)  +  CHp*EHD ; 

%  Filling  in  Negative  k___z  Spectrum  Using  Symmetry 
TEz (nkzp+l:nz)=flipud(TE2 (2: nkzp-1) ) ; 
THp(nk2p+l:nz)=flipud(THp(2:nk2p-D ) ; 

%  Plotting  Magnitudes  of  Transfer  Functions 
%  proper  fftshift  operation 

Ppl  (1  ;nkzp-*2)=abs  (TEz  (nk2p+l:nz) ) ;  Ppl  (nkzp-1  :nz)=abs  (TEz  (l:nkzp)  )  ; 

Pp2  (1  :nkzp-2)=abs  (THp  (nkzp+1  :nz) )  ;  Pp2  (nkzp-1  :nz)=abs  (THp  (1  :nk2p) ) ; 
clf  reset;  plot (kz, Ppl, ’g’ , kz, Pp2, ’r: ’ ) 
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title  ([ 'Transfer  Functions:  |TEz(k_z)|  (solid);'... 

'  iTHp(k_z)|  (dash);  Rhol=' ,num2str (rpl) , . . . 

'm;  Rho2=' ,n\im2str  (rp2) , 'in;  CEz=' ,num2str  (CEz) ,  .  . . 

';  CHp=' ,num2str (CHp) ] ) 

xlabel('k_z  (rad/m)');  ylabel ( 'Magnitude' ) ;  figure(l) 
hep y=input ( 'Enter  1  for  Hard  Copy:  '); 
if  ~iseinpty  (hepy) ,  if  hepy  ==  1,  print;  end;  end 

%  Propagating  Spectra  and  Inverse  Transforming  Fields 
Ezf2=TEz.*Ezf (:,ml) ;  Hpf2=THp. *Hpf ( : ,ml) ;  Erf2=THp.*Erf ( : ,ml) ; 

%  Use  fft  as  inverse  to  ifft  for  k_z  polarity  reversal  in  transforms 
Ez2=fft (Ezf2) /nz;  Hp2=f ft (Hpf 2) /nz;  Er2=fft(Erf2)/nz; 

Pr2=-pi*rp2.*real (E22.*conj (Hp2) ) ;  %  rho-directed  power/dz 

Pz2=+pi*rp2.*real (Er2.*conj (Hp2) ) ;  %  z-directed  power/dz 

Prf2=-rp2 .*real (Ezf2 .*conj (Hpf2) ) /2; %  Power  Spectral  Densities 
Pzf2=+rp2.*real (Erf2.*conj (Hpf2) ) /2; 

%  Computing  RMS  Errors  in  Spatial  Fields  at  Rho2 

DEz=Ez ( : ,m2) -Ez2;  Ezerr=100*sqrt ( (DEz' *DEz) / (Ez ( : ,m2) ' *Ez ( : ,m2) ) ) ; 
DHp=Hp ( : ,m2) -Hp2;  Hperr=100*sqrt ( (DHp' *DHp) / (Hp ( : ,m2) ' *Hp (:,m2))); 
DEr=Er(:,m2)-Er2;  Ererr=100*sqrt ( (DEr' *DEr) / (Er ( : ,m2) '*Er (:,m2) ) ) ; 
DPr=Pr ( : ,m2) -Pr2;  Prerr=100*sqrt ( (DPr ' *DPr) / (Pr ( : ,m2) ' *Pr ( : ,m2) ) ) ; 
DPz=Pz(:,m2)-Pz2;  Pzerr=100*sqrt ( (DPz' *DPz) / (Pz ( : ,m2) '*Pz(:,m2) ) ) ; 

%  Plotting  Spatial  and  Spectral  Fields  at  rho=rp2 
clf  reset 

subplot (2,1,1),  plot (zee, abs (Ez ( : ,m2) ) , ' g' , zee, abs (Ez2) , ' r : ' ) 
title (['E_2  at  Rho=' ,num2str (rp2) , 'm  Propagated  from  Rho=',... 
num2str (rpl) , 'm;  CEz=' ,num2str (CEz) , ' ;  RMS  Error=',... 
num2str (Ezerr) ,'%']) 
xlabelCz  (m)  ' ) ;  ylabel  ( 'V/m' ) 

Ppl (1 :nkzp-2) =abs (Ezf (nkzp+1 :nz,m2) ) ;  Ppl (nkzp- 
1 :nz) =abs (Ezf (1 :nkzp,m2) ) ; 

Pp2  (1 :nkzp-2) =abs (Ezf 2 (nkzp+1 :nz) ) ;  Pp2 (nkzp-1 :nz) =abs (Ezf 2 (1 :nkzp) ) ; 

subplot (2,1,2),  plot (kz, Ppl, ' g' , kz, Pp2, 'r: ' ) 

titleCPlot  Key:  Exact  (solid)  Propagated  (dash)') 

xlabel('k_z  (rad/m)');  ylabel  (' |  FFT  | ')  ;  figured) 

hep y=input ( 'Enter  1  for  Hard  Copy:  '); 

if  ~isempty (hepy) ,  if  hepy  ==  1,  print;  end;  end 

clf  reset 

subplot (2,1,1),  plot (zee, abs (Hp ( : ,m2) ) , 'g' , zee, abs (Hp2) , ' r : ' ) 
title ([ 'H_{phi}  at  Rho=' ,num2str (rp2) , 'm  Propagated  from  Rho=',... 
num2str (rpl) , 'm;  CHp=' ,num2str (CHp) , ' ;  RMS  Error=',... 
num2str (Hperr) ,'%']) 
xlabelCz  (m)  ' )  ;  ylabel  ( 'V/m' ) 

Ppl (1 :nkzp-2) =abs (Hpf (nkzp+1 :nz, m2) ) ;  Ppl (nkzp- 
l:nz)=abs (Hpf (l:nkzp,m2) ) ; 

Pp2 (l:nkzp-2)=abs (Hpf2 (nkzp+l:nz) ) ;  Pp2 (nkzp-1 :nz) =abs (Hpf2 (l:nkzp) ) ; 

subplot (2,1,2),  plot (kz, Ppl, 'g' , kz,Pp2, 'r : ' ) 

titleCPlot  Key:  Exact  (solid)  Propagated  (dash)') 

xlabel('k_z  (rad/m)');  ylabel (' | FFT | ') ;  figure (1) 

hcpy=input ('Enter  1  for  Hard  Copy:  '); 

if  ~isempty (hepy) ,  if  hepy  ==  1,  print;  end;  end 

clf  reset 

subplot (2,1,1),  plot (zee, abs (Er ( : ,m2) ) , 'g' , zee, abs (Er2) , ' r: ' ) 
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title  {[ 'E_{rho}  at  Rho=\niim2str  {rp2) , ’m  Propagated  from  Rho=’,*.. 
num2str  (rpl) , ’m;  CEr=CHp=’  ,nuia2str  (CHp) ,  * ;  RMS  Error=’ ,  . . . 
niim2str  (Ererr)  /’%’]) 
xlabel(*z  (m)  ’ ) ;  ylabel  ( ’ V/m’ ) 

Ppl  ( 1 : nkzp-’2 )  =abs  (Erf  (nkzp+1  :nz,in2) ) ;  Ppl  (nkzp- 
1 :nz) =abs (Erf (1 :nkzp,m2) ) ; 

Pp2 (1 :nkzp-2) =abs (Erf2 (nkzp+1 :nz) ) ;  Pp2 (nkzp-1 :nz) =abs (Erf2 (1 :nkzp) ) ; 

subplot (2, 1,2) ,  plot (kz, Ppl, ’g’ , kz, Pp2,  ^r : ' ) 

title(’Plot  Key:  Exact  (solid)  Propagated  (dash)’) 

xlabel(’k_z  (rad/m)’);  ylabel ( N FFT I* ) ;  figure (1) 

hep y=input (’ Enter  1  for  Hard  Copy:  ’); 

if  -'isempty  (hepy) ,  if  hepy  ==  1,  print;  end;  end 

clf  reset 

subplot (2,1,1),  plot (zee, Pr ( : ,m2) , * g* , zee, Pr2, ’r : * ) 

title ([ ’P_{rho}  at  Rho=’ , num2str (rp2) , ’m  Propagated  from  Rho=',... 

num2str  (rpl) , ’m;  CEz=’  ,num2str  (CEz)  ,  ’ ;  CHp='  ,nxjm2str  (CHp) _ _ 

’;  RMS  Error=’ ,num2str (Prerr) , ’%’ ] ) 
xlabeK’z  (m)’);  ylabel  ( ’W/m’ ) 

Ppl  (l:nkzp-2)=Prf (nkzp+1 :nz, m2) ;  Ppl (nkzp“l:nz)=Prf (1 :nkzp,m2) ; 

Pp2 (l:nkzp-2)=Prf2 (nkzp+1 :nz) ;  Pp2 (nkzp-l:nz)=Prf2 (l:nkzp) ; 
subplot (2,1,2),  plot (kz, Ppl, * g’ , kz, Pp2,  ’r : ’ ) 
title(’Plot  Key:  Exact  (solid)  Propagated  (dash)’) 

xlabel(’k_z  (rad/m)’);  ylabel ( ’P_{rho}  PSD’);  figure (1) 
hcpy=input ( ’Enter  1  for  Hard  Copy:  ’); 
if  -isempty (hepy) ,  if  hepy  ==  1,  print;  end;  end 

clf  reset 

subplot (2,1,1),  plot (zee, Pz ( : ,m2) , ’g’ , zee, Pz2, ’r : ’ ) 
title  ([’P_z  at  Rho=’ ,niim2str  (rp2) , ’m  Propagated  from  Rho=’,... 
num2str (rpl) , ’m;  CEr=CHp=’ , num2str (CHp) , . . . 

’;  RMS  Error=’ ,num2str (Pzerr) , ’%’ ] ) 
xlabeK’z  (m)  ’ ) ;  ylabel  ( ’ W/m’ ) 

Ppl (1 :nkzp-2) =Pzf (nkzp+l:nz,m2) ;  Ppl (nkzp-1 :nz) =Pzf (l:nkzp,m2) ; 

Pp2 (l:nkzp-2)=Pzf2 (nkzp+l:nz) ;  Pp2 (nkzp-1 :nz)=Pzf2 (l:nkzp) ; 

subplot (2,1,2),  plot (kz, Ppl, ’g’ , kz, Pp2, ’ r : ’ ) 

title (’Plot  Key:  Exact  (solid)  Propagated  (dash)’) 

xlabel(’k_z  (rad/m)’);  ylabel (’P_z  PSD’);  figure (1) 

hcpy=input (’ Enter  1  for  Hard  Copy:  ’); 

if  -'isempty  (hepy) ,  if  hepy  ==  1,  print;  end;  end 

end;  end 

if  nplot==9, 

if  isempty(rl),  error(’Enter  Rhol  First’);  end 

r2=input ( ’Enter  Rho2  for  Endpoint  of  Propagation  Error  Calculation:  ’); 

if  r2==rl,  error (’Enter  Rho2  Different  Then  Rhol’);  end 

%  Selecting  closest  grid  radius 

[rm,  m2]=min(abs  (rho-r2)  )  ; 

if  m2  >  ml,  M=ml:m2;  else,  M=m2:ml;  end 

rp=rho (M) ;  nrp=length (rp) ; 

FE2=zeros (nkzp,nrp) ;  FH2=FE2;  TEz=zeros (nz,nrp) ;  THp=TEz; 

%  Forming  denominator  of  radial  transfer  functions 
for  m=l:nrp; 

[H,  DH]=dhankel (0, krsup*rp (m) ) ; 

FE2 (l:nkzsup,m)=H;  FH2 (1 :nkzsup,m) =DH; 
if  --isempty  (krsub) ,  %  subluminal  spectra 

[K,  DK]=dmodbes  (0,  krsub*rp  (m) )  ; 
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FE2 (nkzsup+1 : nkzp, m) =K; 

FH2  (nkzsup+1: nkzp,  la)  =DK; 
end;  end 

FEE=FEl*ones {l,nrp) ;  FHH=FHl*ones {l,nrp) ;  Ppl=zeros (nz,nrp) ;  Pp2=Ppl; 
while  2 

CEz=input ( 'Enter  Optimal  Estimator  Constant  for  E_z  (CEz  <  0  to  exit): 
'); 

if  CEz  <  0,  break;  end 

CHp=input (' Enter  Optimal  Estimator  Constant  for  H_phi  &  E_rho:  '); 

%  Computing  Riad's  Optimal  Deconvolution  Propagation  Transfer  Functions 
TEz(l:nkzp, :)=FE2.*conj (FEE) ./{FEE.*conj (FEE)  +  CEz*EEl) ; 

THp(l:nkzp, :)=FH2.*conj (FHH) ./(FHH.*conj (FHH)  +  CHp*EHl) ; 

%  Filling  in  Negative  k_z  Spectrum  Using  Symmetry 
TEz (nkzp+l:nz, : )=flipud(TEz (2:nkzp-l, : ) ) ; 

THp (nkzp+1 :nz, : )=flipud(THp(2:nkzp-l, : ) ) ; 

%  Plotting  Magnitudes  of  Transfer  Functions 

%  proper  fftshift  operation 

Ppl (1 :nkzp-2, : ) =abs (TEz (nkzp+1 :nz, : ) ) ; 

Ppl (nkzp-l:nz, : ) =abs (TEz (l:nkzp, : ) ) ; 

Pp2  U:nkzp-2, : )=abs (THp (nkzp+1 :nz, : ) ) ; 

Pp2 (nkzp-l:nz, : )=abs (THp(l:nkzp, : ) ) ; 

clf  reset;  mesh (rp, kz,  Ppl) 

title ( [ ' Deconvolution  Function:  1 TEz (kz, rho) |  with  ' . . . 

'CEz=' ,num2str (CEz) ,  ' ;  Propagate  from  Rhol=' ,num2str (rpl) , 'm' ] ) 
ylabel('k_z  (rad/m)');  xlabel('rho  (m) ' ) ;  view (135, 40) ;  figure (1) 
hcpy=input (' Enter  1  for  Hard  Copy:  '); 
if  ~isempty (hcpy) ,  if  hcpy  ==  1,  print;  end;  end 

clf  reset;  mesh(rp,  kz,  Pp2) 

title ([' Deconvolution  Transfer  Function:  | THp (kz, rho) I  for  ' — 

'CHp=' ,num2str (CHp) ,  ' ;  Propagate  from  Rhol=' ,num2str (rpl) , 'm' ] ) 
ylabel('k_z  (rad/m)');  zlabeK'rho  (m)  ' )  ;  view(135, 40)  ;  figured) 
hcpy=input ( 'Enter  1  for  Hard  Copy:  '); 
if  ~isempty (hcpy) ,  if  hcpy  ==  1,  print;  end;  end 

%  Propagating  Spectra  and  Inverse  Transforming  Fields 
Ezf2=TEz.* (Ezf ( : ,ml) *ones (l,nrp) ) ; 

Hpf2=THp.* (Hpf ( : ,ml) *ones (l,nrp) ) ; 

Erf2=THp.* (Erf ( : ,ml) *ones (l,nrp) ) ; 

Ez2=fft(Ezf2)/nz;  Hp2=f ft (Hpf2) /nz;  Er2=fft(Erf2)/nz; 

Pr2=-pi* (ones (nz, 1) *rp) .*real (Ez2 .*conj (Hp2) ) ;  %  rho-directed  power/dz 

Pz2=+pi* (ones (nz, 1) *rp) •*real (Er2 .*conj (Hp2) ) ;  %  z-directed  power/dz 

%  Computing  RMS  Errors  in  Spatial  Fields  at  Rho2 
DEz=Ez ( : , M) -Ez2 ;  DHp=Hp ( : , M) -Hp2 ;  DEr=Er ( : , M)  -Er2 ; 

DPr=Pr ( : , M) -Pr2 ;  DPz=Pz ( : , M) -Pz2 ; 
for  m=l:nrp; 
m2=M (m) ; 

Ezerr(m)=100*sqrt( (DEz(:,m) '*DEz(:,m) )/ (Ez(:,m2) '*Ez(:,m2) ) ) ; 
Hperr(m)=100*sqrt( (DHp(:,m) ' *DHp ( : ,m) ) / (Hp ( : ,m2) '*Hp(:,m2) ) ) ; 

Ererr (m)=100*sqrt ( (DEr ( : ,m) '*DEr ( : ,m) ) / (Er ( : ,m2) '*Er ( : ,m2) ) ) ; 

Prerr(m)  =100*sqrt{  (DPr{:,m)  '*DPr(:,m) )/  (Prd, m2)  '*Pr(:,m2) ) )  ; 
Pzerr(m)=100*sqrt ( (DPz(:,m) ' *DPz ( : ,m) ) / (Pz ( : ,m2) '*Pz(:,m2) ) ) ; 
end 
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clf  reset;  plot (rp,E2err, ’ g' , rp,Hperr, ’r: ’ , rp,Ererr/ * c — ’ ) 
xlabel{’Rho  (m)  ’ ) ;  ylabel(^RMS  Error  (%)’)/ 
v=  axis;  v(3)=0;  axis  (v)  ;  figured) 

title ([ ’Error  in  E_z  (solid),  H_{phi}  (short  dash)  and  E_{rho}’. 
’  (long  dash)  Propagated  from  Rhol=’ ,n\im2str  (rpl) ,  . . . 

’ ;CEz=’ ,num2str (CEz) , ’ ;  CHp=’ ,num2str (CHp) ] ) 
hep y=input ( ’Enter  1  for  Hard  Copy:  ’); 
if  -'isempty  (hepy) ,  if  hepy  ==  1,  print;  end;  end 

clf  reset;  plot (rp, Prerr, ’ g’ , rp, Pzerr, ’r; ’ ) 
xlabel(’Rho  (m)  ’ ) ;  ylabel(’RMS  Error  (%)’) 
v=  axis;  v(3)=0;  axis (v) ;  figure (1) 

title  ([ ’Error  in  P_r  (solid)  and  P__2  (dash)  Propagated  from  Rhol 
num2str (rpl) , ’ ;CEz=’ , num2str (CEz) , ’ ;  CHp=’ ,num2str (CHp)  ] ) 
hcpy=input ( ’Enter  1  for  Hard  Copy:  ’); 
if  is  empty  (hepy) ,  if  hepy  ==  1,  print;  end;  end 

end;  end;  end 
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%  EHFLDTST9.M  tests  dRhoEHFld  function  routine  by  comparing  computation 
%  for  sinusoidal  I(z)=Im  sin  k{h-l2|)  at  selected  cylindrical 
%  radius*  (Dipole  at  the  origin  with  sinusoidal  Current  Distribution) . 

%  Using  exact  integration  results  for  this  assumed  current 
%  given  in  Jordan  &  Balmain  pp  333-337  and  Elliot  pp  329-332  for 
%  the  Exact  Fields  at  radius  1  and  radius  2  (>  radius  1)  and  the 
%  Computed  Ez  and  Hp  fields  at  radius  2  given  the  fields  at  radius  1. 

%  Adapted  from  EHFLDTST.M  by  M.A.  Morgan  6/6/97 
%  by  D.G.  Steenman  12/28/98 

%  Mod  9  of  EHFLDTST:  Calls  dRhoEHFld4  which  uses  different  field 
%  calculations. 

%  Calls  dRhoEHFld4  (which  calls  intphi3  and  phisquare) 
clear  all 

L=input ( ’Enter  dipole  length  L  (meters):  ’);  h=L/2; 

f=input ( ’Enter  frequency  (MHz) :  ’ ) ; 

wl=300/f;  k=2*pi/wl;  kh=k*h;  Im=l;  ckh=cos(kh); 

Invk  =  1/k; 

disp([’k*h  =  ’ ,num2str (kh) ] ) ; 
disp([’l/k  =  ’ ,num2str (Invk) ] ) ; 

while  1 


zmax=input  ( ’Enter  zmax  (m)  for  z=[-zmax  .*.  zmax]  plots  (0  to  stop):  ’); 
if  -'isempty  (zmax) ,  if  zmax  =-  0,  break;  end;  end 
M=input ( ’Enter  #  z-points  for  [-zmax  . . *  zmax] :  ’ ) ; 
dz=2*zmax/ (M-1) ;  2= (-zmax: dz : zmax) ; 

rhol=input ( ’Enter  rhol  (m)  for  plot  (0  to  stop):  ’); 
if  '-isempty  (rhol) ,  if  rhol  ===  0,  break;  end;  end 

rho2min=input  ( ’Enter  rho2  min  (m)  for  [rho2min  ...  rho2max]  (0  to  stop): 
;)  ; 

if  '^isempty  (rho2min)  ,  if  rho2min  =-  0,  break;  end;  end 

rho2max=input  ( ’Enter  rho2  max  (m)  for  [rho2min  ...  rho2max]  (0  to  stop): 
’); 

if  ^isempty  (rho2max)  ,  if  rho2max  ==  0,  break;  end;  end 

Nrho  =  input  (’  Enter  #  rho2-points  for  [rho2min  ...  rho2max]  ;  ’); 

dr=  (rho2max-rho2min)  /  (Nrho-1) ;  DRHO=(rho2min:dr:rho2max) ; 

NpCalc  =  ceil (pi*2*rhol/dz) 

Nplow=input ( ’Enter  Min  number  of  Phi  segments  around  a  thin  ring  (radius 
1)  :  M; 

if  NpCalc>Nplow, 

Np=NpCalc 

else, 

Np~Nplow 

end 

%  Exact  Field  Calculation  at  Radius  1 
2l=z-h;  z2=2+h; 

r=sqrt (rhol*rhol+z.*2) ;  E=exp (-j*k*r) ; 

Rl=sqrt (rhol*rhol+zl.*2l) ;  El=exp(-j*k*Rl) ; 

R2=sqrt (rhol*rhol+z2.*22) ;  E2=exp {-j*k*R2) ; 
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Ezl=-j*30*lm* (El./Rl+E2./R2-2*ckh*E./r)  ; 

Hpl=(j*Im/ (4*pi) ) * (El+E2-2*ckh*E) /rhol; 

Jz=Hpl;  Mp=Ezl; 

%  Set  Up  Loop  to  do  sets  of  rho2  Field  Calculations 
EZRH02  =  zeros (length (z) , length (DRHO) );  ■ 

HPRH02  =  zeros (length (z) , length (DRHO) )  ; 

EZRH02N  =  zeros (length (z) , length (DRHO) )  ; 

HPRH02N  =  zeros (length (z) /length (DRHO) )  ; 

EZRH02NT  =  zeros (length (z) , length (DRHO) )  ; 

HPRH02NT  =  zeros (length (z) /length (DRHO) )  ; 

for  a  =  1 : (length (DRHO) ) / 
rho2  =  DRHO (a); 

%  Exact  Field  Calculation  at  Radius  2 
rrho2=sqrt (rho2*rho2+z . *z) ;  Erho2=exp (-j*k*rrho2) ; 

Rlrho2=sqrt (rho2*rho2+zl.*zl)  ;  Elrho2=exp(-j*k*Rlrho2) ; 

R2rho2=sqrt (rho2*rho2+z2.*z2) ;  E2rho2=exp (-j*k*R2rho2) ; 

Ez2=- j *30*Iia*  (Elrho2 . /Rlrho2+E2rho2 ,  /R2rho2-2*ckh*Erho2 .  /rrho2) ; 
Hp2=(j*Im/ (4*pi) )* (Elrho2+E2rho2-2*ckh*Erho2)/rho2; 

EZRH02(:/a)=Ez2.'; 

HPRH02 ( : / a) =Hp2 . ' ; 

%  Numerical  Integration  using  dRhoEHFld  function  for  Radius  2 
[Ez2N,  Ez2NT/  Hp2N/  Hp2NT]=dREHF4(Np/k/rhol/rho2/Z/Z/Z/Z/l/l/JZ/Mp); 
EZRH02N(:/a)=Ez2N. ' ; 

HPRH02N ( : / a) =Hp2N. ' ; 

EZRH02NT (: /a)=Ez2NT. ' ; 

HPRH02NT ( : / a) =Hp2NT . ' ; 

%clf  reset; 
figure (2*a-l) 

plot  (Z  /  real  (Ez2)  /  '  q'  ,z,  real  (Ez2N) '.r'  ,z,  real  (Ez2NT)  /  '  -  .y'  /  . . . 

Z/  imag  (Ez2)  ,  '  q'  ,z,  imag  (Ez2N) -.r'  ,z,  imag  (Ez2NT) ,  '  -  .y' ,  .  . . 
z,  abs (Ez2) / 'q' ,z, abs (Ez2N) , ' : r ' / Z/ abs (Ez2NT) / '-.y' ) 
title (['E_z:  Exact  (solid);  Integ  (dot);  NullTest  (dashdot)  at 
{\rho}2='  / num2str (rho2) / . - . 

'  m  from  {\rho}l=  ’ /num2str (rhol)  /  '  m  for  L=' , num2str (L) / ’  m; 

N_z= ' / num2str (M) / . . . 

';  N_{\phi}=' / int2str  (Np)  / ' ;  f=’ /num2str (f )  /  '  MHz']); 
xlabel(['z  (m)  /  Density  N_z  =  ' /num2str  (M/ (2*zmax) )  ] ) ;  ylabelCReal  and 
Imag  E_z  (V/m) ' ) ; 
grid;  orient  landscape; 

%hcpy=input (' Print  Hard  Copy  ?  (Y/N) :  '/'s'); 

%if  ~isempty (hcpy) /  if  hcpy=='Y' |hcpy=='y' /Orient  landscape;  print;  end; 
end 

%clf  reset; 
figure (2*a) 

plot (Z  /  real (Hp2) ,  'q' ,z,  real  (Hp2N) / ' :r ' / Z/ real (Hp2NT)  /  '- .y' / . . . 

Z/  imag (Hp2) / ' q' ,z, imag  (Hp2N) \r' ,z, imag (Hp2NT) ,  ' - .y' ,  . . . 
z,  abs (Hp2) / 'g' / Z/ abs (Hp2N) w' ,z,  abs (Hp2NT) ,  '-.y' ) 
title ([ 'H_{phi} :  Exact  (solid);  Integ  (dot);  NullTest  (dashdot)  at 
{\rho}2='  / num2str (rho2) ,  . . . 
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'  m  from  {\rho}l=  ' , num2str  (rhol) , '  m  for  L=' ,niim2str  (L) ,  '  m; 

N_z= ’ / num2str (M) / . . . 

N_{\phi}=*,int2str{Np), f=’ ,num2str (f) , '  MHz']); 
xlabel{['z  (m) ,  Density  N_z  =  '  ,num2str  (M/ (2*zmax) )  ] ) ;  ylabeK'Real  and 
Imag  H_{phi}  (A/m)*); 
grid;  orient  landscape; 

%hcpy=input ( ' Print  Hard  Copy  ?  { Y/N) ;  ' , ' s ' ) ; 

%if  ~isempty(hcpy) ,  if  hcpy=='Y' |hcpy=='y' ,  orient  landscape;  print; 
end;  end 

end 

d  =  2*length(DRH0); 

%d=0; 

[NRH02G,ZEEG]=meshgrid{DRH0./rhol,  z) ; 

RErMagH  =  { {abs (HPRH02-HPRH02N) ) ) ./ { (abs (HPRH02) ) ) ; 

RErMagE  =  ( (abs (EZRH02-EZRH02N) ) ) . / ( (abs (EZRH02) ) ) ; 

LSErMagH  =  sqrt ( (sum( (abs (HPRH02- 

HPRH02N) )  .''2)  )  ./  (sum(  (abs  (HPRH02) )  .''2) ) )  *100; 

LSErMagE  =  sqrt ( (sum ( (abs (EZRH02- 

EZRH02N) )  .''2) )  ./  (sum(  (abs  (EZRH02) )  .^2)  ) )  *100; 

RErMagHNT  =  ( (abs (HPRH02-HPRH02NT) ) ) ./ ( (abs (HPRH02) ) ) ; 

RErMagENT  =  ( (abs (EZRH02-EZRH02NT) ) ) . / ( (abs (EZRH02) ) ) ; 

LSErMagHNT  =  sqrt ( (sum( (abs (HPRH02- 
HPRH02NT) )  .^2) )  ./  (sum(  (abs(HPRH02) )  .''2)))  *100; 

LSErMagENT  =  sqrt (  (sum(  (abs (EZRH02- 
EZRH02NT) )  .'‘2))  ./  (sum(  (abs  (EZRH02) )  .''2) ) )  *100; 


figure (d+1) 

surf (NRH02G, ZEEG, RErMagH) 

title(['  Relative  Error  in  Exact  to  Num  Soln  for  Mag(H_{\phi}  ({\rho} 
2));  {\rho}  1=  ' ,num2str (rhol) ,  '  (m)  ']) 
xlabeKC  {\rho}  2/{\rho}  1  ;  N_{\rho}=  '  ,num2str  (Np)  ] )  ; 
ylabel(['  z  (m)  ;  N_z  =  ' ,num2str (M) ,  '  Density  N_z  = 

'  ,num2str  (M/  (2*zmax) )  ] ) ; 
orient  landscape; rotateSd; 

figure (d+2) 

surf (NRH02G, ZEEG, RErMagHNT) 

title(['  Relative  Error  in  Exact  to  Num  Null  Test  Soln  for  Mag(H_{\phi} 
({\rho}  2));  {\rho}  1  =  ' , num2str (rhol) , '  (m)  ']) 

xlabeKC  {\rho}  2/{\rho}  1  ;  N_{\rho}=  '  ,num2str  (Np)  ] )  ; 
ylabeKE'  z  (m)  ;  N_z  =  ' , num2str  (M) ,  '  Density  N_z  = 

'  ,num2str  (M/  (2*zmax) )  ] ) ; 
orient  landscape ;rotate3d; 

figure (d+5) 

surf (NRH02G, ZEEG, RErMagE) 

title (['  Relative  Error  in  Exact  to  Num  Soln  for  Mag(E_z  ({\rho}  2) ) ; 
{\rho}  1  =  '  ,nina2str  (rhol) ,  '  (m)  ']) 

xlabel(['  {\rho}  2/{\rho}  1  ;  N_{\rho}=  ' ,num2str (Np) ] ) ; 
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ylabel  ( [ '  z  (m)  ;  N_z  =  '  ,niim2str  (M) ,  '  Density  N_z  = 

’,n™2str(M/  (2*zmax)  )  ]  )  ; 
orient  landscape;  rotateSd; 

figure (d+6) 

surf (NRH02G, ZEEG,RErMagENT) 

title(['  Relative  Error  in  Exact  to  Num  Null  Test  Soln  for  Mag(E_z 
({\rho}  2));  {\rho}  1  =  ’ ,num2str (rhol) , '  (m)  ']) 

xlabel ( [ '  {\rho}  2/{\rho}  1  ;  N_{\rho}=  ' ,num2str (Np) ] ) ; 
ylabel  (['  z  (m)  ;  N_z  =  '  ,nuia2str  (M) ,  ’  Density  N_z  = 

' ,  nuirL2str  (M/  (2*zmax)  )  ] )  ; 
orient  landscape;  rotateSd; 

figure (d+7) 

plot (DRHO, LSErMagH, ’ r ' , DRHO, LSErMagE, ' — 
g',DRHO,LSErMagHNT, ' :y’ ,DRHO,LSErMagENT, '-.y') 

title(['  Least  Squares  Error  for  Mag(H_{\phi})  and  Mag(E_z)  ;  {\rho}  1 
' ,num2str (rhol) , '  m  ']) 

xlabel ([ *  {\rho}  2  ;  N_{\rho}=  ' ,num2str (Np) , '  Density  N_z  = 

' ,num2str (M/ (2*zmax) ) ] ) ; 

ylabel ( '  Least  Squares  Error  (%);  MagH  (solid),  MagE  (dash),  MagH 
NullTest  (dots),  MagE  NullTest  (dash-dot)  '); 
grid;  orient  landscape; 


end 
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function  [EzJ,  EzM/  HpJ,  HpM,  Tezj,  Tezm,  Thpj,  Thpia,  zOffset]  = 
dREH JMA (Np,  k, rho 1 , rho2 , z 1 , z2 , Jz , Mp ) 

%  Method  A  for  Niimerically  Calculating  the  Fields  at  Rho2. 

%  Calculate  the  EzJ,  EzM,  HpJ,  and  HpM  fields  on  a  cylindrical  surface 
of 

%  radius  2  due  to  the  equivilent  Electric  and  Magnetic  Currents 
%  on  a  cylindrical  surface  of  radius  1.  (Radius  2  >  Radius  1), 

%  Fields  are  calculated  by  integrating  a  constant  magnitude  current 
%  around  a  thin  ring  of  radius  1  and  summing  all  the  contributing 
%  pieces  of  fields  from  each  element  in  z. 

%  The  Transforms  are  an  interim  step  to  finding  the  fields. 

%  Includes  ability  to  calculate  fields  on  a  z2  smaller  than  zl. 

%  Adapted  from  M.A.  Morgan  notes  10  Dec  98 

%  Calls  intphiS  and  phisquare2 

%  Calculate  the  Phi  Integrations  for  Nphi,  zd,  Rhol,  and  Rho2 
zd  =  zl-zl (1) ; 

[u2,u3,u4,u5,  v2,  v3,  v4,  v5,dv3,  dv4,  dv5,  w3,  w4,  w5]  = 
intphi3  (Np,  k,  rhol,  rho2,  zd) ; 

%  Set  Up  Integration  Matrices  from  Calculated  Nphi  segmented  phi 

%  integrations  around  a  thin  ring  of  radius  rhol 

znl=length (zl) ; 

zn2=length(z2) ; 

izmax  =  find (zl> (max (z2) )) ; 

zOffset=length(zl (izmax) ) ; 

U2=  phisquare2 (u2, znl, zn2, zOf fset) ; 

U3=  phisquare2 (u3, znl, zn2, zOf fset) ; 

U4=  phisquare2 (u4, znl, zn2, zOffset) ; 

U5=  phisquare2 (u5, znl, zn2, zOffset) ; 

V2=  phisquare2 (v2, znl, zn2, zOffset) ; 

V3=  phisquare2 (v3, znl, zn2, zOffset) ; 

V4=  phisquare2 (v4, znl, zn2, zOffset) ; 

V5=  phisquare2 (v5, znl, zn2, zOffset) ; 

DV3=  phisquare2 (dv3,  znl, zn2, zOffset) ; 

DV4=  phisquare2 (dv4, znl, zn2, zOffset) ; 

DV5=  phisquare2 (dv5, znl, zn2, zOffset) ; 

W3=  phisquare2 (w3, znl, zn2, zOffset) ; 

W4=  phisquare2 (w4, znl, zn2, zOffset) ; 

W5=  phisquare2 (w5, znl, zn2, zOffset) ; 


%  Calculate  the  Components  due  to  the  Different  Current  Components 
%  Method  A  for  solving  the  fields 

dz  =  abs (zl (2) ~zl (1) ) ; 

%  Hp  due  to  Jz  (Hp  at  rho2  due  to  Hp  at  rhol) 

hpj  =  (rho2*U3  -  rhol*V3)  +  j*k*(rho2*U2  -  rhol*V2); 

Thpj  =  (rhol*dz/ (4*pi) ) *hpj ; 

Hpj  =  Jz*Thpj; 
clear  hpj ; 

%  Ez  due  to  Jz  (Ez  at  rho2  due  to  Hp  at  rhol) 
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Rezj  =  {{2+k*k*rho2*rho2)*U3  -  3*rho2*rho2*U5  - 
( (rhol/rho2)+2*k*k*rhol*rho2)*V3  • • . 

+  6*rhol*rho2*V5  +  k*k*rhol*rhol*W3  -  3*rhol*rhol*W5) ; 

Imezj  =  j*k*(2*U2  -  3*rho2*rho2*U4  - (rhol/rho2) *V2  +  6*rhol*rho2*V4  - 
3*rhol*rhol*W4) ; 

Tezj  =  (rhol*dz*30/ ( j*k) ) * (Rezj  +  Imezj); 

EzJ  =  Jz*Tezj; 
clear  Rezj  Imezj; 

%  Hp  due  to  Mp  (Hp  at  rho2  due  to  Ez  at  rhol) 

Rhpiti  =  (3*rhol*rho2*U5  -  k*k*rhol*rho2*U3 
+ (2+k*k*rhol*rhol+k*k*rho2*rho2) *V3  ... 

+  k*k*DV3  -  3* (rhol*rhol+rho2*rho2) *V5  -3*DV5  -k*k*rhol*rho2*W3  + 
3*rhol*rho2*W5) ; 

Imhpm  =  j*k* (3*rhol*rho2*U4  +  2*V2  -  3* (rhol*rhol+rho2*rho2) *V4  -  3*DV4 
+  3*rhol*rho2*W4) ; 

Thpm  =  (rhol*dz/  ( j*k*4*120*pi*pi) )  *  {Rhpm+ Imhpm)  ; 

HpM  =  Mp*Thpm; 
clear  Rhpm  Imhpm; 

%  Ez  due  to  Mp  (Ez  at  rho2  due  to  Ez  at  rhol) 

ezm  =  (rho2*V3  -  rhol*U3)  +  j*k*(rho2*V2  -  rhol*U2) ; 

Tezm  =  (rhol*dz/ (4*pi) ) *ezm; 

EzM  =  Mp*Tezm; 
clear  ezm; 
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function  [U1,U2,U3,U4,U5,DU3,DU4,DU5, V1,V2,V3,W3,W4,W5]  = 
intphi4 (N, k, rhol, rho2, z2) 

%  Program  uses  trapazoidal  rule  to  evaluate  integrals  on  phi=0  to  2pi 
circular  contour 

%  of  radius  rhol  in  the  x-y  plane  with  field  point  at  (rho2,  z2,phi=0)  of 
even  functions: 

%  un(phi)=  (1/R^n)  exp(“jkR) 

%  vn(phi)=  cos (phi)  un(phi) 

%  Dun (phi) =  (z2^2)  un(phi) 

%  and  wn(phi)=  cos"^2(phi)  un(phi) 

%  where  R=sqrt  (rhol^2+rho2'^2-2rhol*rho2*cos  (phi) +z2'"2) 

%  Adapted  from  IntPhi.m  and  IntPhi2.m  by  M.A*  Morgan  12/28/98 

dp=pi/(N-l);  phi=0:dp:pi;  %  double  the  integral  from  0  to  pi 
a=rhol*rhol+rho2*rho2+z2.^2;  b=2*rhol*rho2;  cp=cos(phi)  ; 

A=a. '  *ones  (1,  length  (cp) ) ;  D=  (z2  .'^2)  .  *  *ones  (1,  length  (cp) ) ; 

CP=ones (length (22) , 1) *cp; 

R=sqrt (A-b*CP) ; 
ul=exp (^j*k*R) ./ (R) ; 

u2=  ul,/R;  u3=  u2,/R;  u4=  u3./R;  u5=  u4./R; 
vl=CP^*ul;  v2=CP**u2;  v3=CP.*u3;  v4=CP**u4;  v5=CP.*u5; 
w3=CP.*v3;  w4=CP.*v4;  w5=CP.*v5; 
du3=D.*u3;  du4=D**u4;  du5=D.*u5; 
c=2*dp;  ds=c*ones (N, 1) ;  ds(l)=c/2;  ds(N)=c/2; 

Ul=ul*ds;  U2=u2*ds;  U3=u3*ds;  U4=u4*ds;  U5=u5*ds; 

Vl=vl*ds;  V2=v2*ds;  V3=v3*ds; 

W3=w3*ds;  W4=w4*ds;  W5=w5*ds; 

DU3=du3*ds;  DU4=du4*ds;  DU5=du5*ds; 
clear  A  D  CP  R  u*  v*  w’^  dv*  ds  cp  phi; 
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function  A  =  phisquare4  (Vecn,Nzl,N22,MC,N,  zOffsetMax,  zOf fsetMin) 

% 

%  Mod4:  (length  of  Intdzl)  *N==  (length  of  d22)  :  Also,  use  MATIAB 
%  ’toeplitz’'  function  instead  of  spdiags, 

% 

%  Mod3:  Takes  different  dzl  |  dz2  sizes  AND 
%  different  zlMax  |  z2Max  values 

%  into  account 

% 

%  PhisquareS  takes  a  vector  un  and  makes  a  rectangular  diagonal  matrix 
%  out  of  its  elements.  ’Vecnd)  ’  is  the  main  diagonal,  'Veen (2)  ’ 

%  is  the  first  diagonal  off  the  main,  both  above  and  below 
%  the  diagonal.  'Veen'  is  a  column  vector.  After  making  a 
%  large  rectangular  matrix,  specific  rows  are  choosen  to  account 
%  for  the  difference  in  dzl  compared  to  dz2. 

%  If  they  are  the  dz's  are  equal,  this  program  makes  a  banded 
%  Toeplitz  matrix. 

%  by  D.G.  Steenman  12/30/98 


%  Set  up  elements  for  each  row 

U1  =  Veen . ' ; 

ANz  =  length  (Veen) ; 

U  =  toeplitz (U1,U1) ; 

%  First,  Deal  with  Offset 

%  Second,  Choose  Selected  Columns  that  correspond  to  zR2  points 
%  Third,  reshape  and  add  N*MC  rows  to  average  the  integrations 
%  Fourth,  squeeze  and  transpose  to  correct  orientation  for  A 

if  (zOffsetMax<=0) , 

U=U(:,  (~zOf fsetMax+1)  :  (ANz+zOffsetMin) ) ; 
else, 

U=U  ( (zOffsetMax+1)  :  (ANz~zOffsetMin) ,  : )  ; 

end 


U==U(:,  (  (0:Nz2-‘l)*N  +  1)); 
U=U. ' ; 

U=reshape(U,Nz2,N*MC,Nzl)  ; 
U=sum(U,2)  ; 

Unsqueeze  (U) ; 

A=U. 

clear  U*  ANz; 
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